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ABSTRACT 

The Solar Dynamic Power Module being developed for Space 
Station Freedom uses a eutectic mixture of LiF-CaF 2 phase 
change material (PCM) contained in toroidal canisters for 
thermal energy storage. Presented herein are the results 
from heat transfer analyses of a PCM containment canister. 
One- and two-dimensional finite-difference computer models 
are developed to analyze heat transfer in the canister 
walls, PCM, void, and heat engine working fluid coolant. 

The modes of heat transfer considered include conduction in 
canister walls and solid PCM, conduction and pseudo - free 
convection in liquid PCM, conduction and radiation across 
PCM vapor filled void regions and forced convection in the 
heat engine working fluid. Void shape, location, growth or 
shrinkage (due to density difference between the solid and 
liquid PCM phases) are prescribed based on engineering 
judgement. The PCM phase change process is analyzed using 
the enthalpy method. The discussion of results focuses on 
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how canister thermal performance is affected by free 
convection in the liquid PCM and void heat transfer. 
Characterizing these effects is important for interpreting 
the relationship between ground-based canister performance 
(in 1-g) and expected on-orbit performance (in micro-g). 

Void regions accentuate canister hot spots and temperature 
gradients due to their large thermal resistance. Free 
convection reduces the extent of PCM superheating and lowers 
canister temperatures during a portion of the PCM thermal 
charge period. Surprisingly small differences in canister 
thermal performance result from operation on the ground and 
operation on-orbit. This lack of a strong gravity 
dependency is attributed to the large contribution of 
container walls in overall canister energy redistribution by 
conduction . 
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NOMENCLATURE 


c = Specific Heat, J/g-K 

A = Void Surface Element Area, cm 2 or constant 

B = Constant 

C5 = Constant 

D = Cooling Fluid Tube Inner Diameter, cm 
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iv = Grid Element Which Contains X v 

iv' = Combined Grid Element 

k = Thermal Conductivity, W/cm-K 

L = Annular Canister Length, cm 

L' = Slab Canister Thickness, cm 

LH = Liquid PCM Vertical Layer Height, cm 

LiF- = Lithium Fluoride-Calcium Fluoride 
CaF 2 

m = He/Xe Mass Flow Rate, g/sec 
MF = Mass Fraction PCM in Element iv 
MFL = PCM Mass Fraction Liquid 

NRS = Number of Void Radiating Surface Elements 
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= Nusselt Number 
= Phase Change Material 
= Prandtl Number 
= Heat flux, W/cm 2 
= Thermal Power, W 
= Radial Coordinate, cm 
= Thermal Resistance, cm 2 -K/W 
= Rayleigh Number 
= Reynolds number 
= Conduction path length, cm 
= Stefan Number 
= Time, sec 
= Temperature, K 
= Thermal Energy Storage 

= Velocity of Void-PCM Interface, cm/sec 
= Overall Heat Transfer Coefficient, W/cm -K 
= Void Volume Fraction 
= Liquid PCM Width, cm 
= One-Dimensional Coordinate, cm 
= Melt or Void Front Position, cm 
= Mushy Zone Liquid PCM Mass Fraction 
= Mushy Zone Liquid PCM Volume Fraction 
= Axial Coordinate, cm 
= Dot Product 
= Gradient Operator 
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Greek Symbols 

a = Thermal Diffusivity, cm 2 /sec 

p = Volumetric Thermal Expansion Coefficient, 1/K 

6 = Canister Wall Thickness, cm 

A = Kronecker Delta Function 

Ar = Radial Grid Spacing, cm 

At = Time Step, sec 

AT = Temperature Difference, K 

Ax = One-Dimensional Grid Spacing, cm 

Az = Axial Grid Spacing, cm 

€ = Emittance 

n = Pi Constant, 2*sin _1 (l) 

p = Density, g/cm 3 

a = Stef an-Boltzmann Constant, 5 . 67051xl0 -12 W/cm 2 -K 4 

v = Kinematic Viscosity, cm 2 /sec 

i|i = Dimensionless Function of Ar and r 

Subscripts and Superscripts 
E = Enhanced 

EFF = Effective 

f = He-Xe Fluid 

i = Inner Radius or Grid Element Index 

j = Void Surface Element 

k = Void Surface Element 

L = Liquid PCM 

m = PCM Melt 

n = Current Time Step 

n+1 = Future Time Step 
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O = Outer Radius 

PCM = Phase Change Material 
rad = Radiation 
S = Solid PCM 
v = Void 
VAP = Vapor 
w = Canister Wall 
WS = Side Wall 
+ = Slightly Above 

= Slightly Below 
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CHAPTER I 


SUMMARY 


Phase change thermal energy storage is a particularly 
attractive approach to meet energy storage requirements for 
the space station Freedom electrical power system. In this 
application, the ability to produce continuous electrical 
power with the intermittent solar source of low earth orbit 
is crucial. The solar dynamic power module proposed for use 
on Freedom incorporates a solid-to-liquid phase change 
material (PCM) encapsulated in multiple, annular containment 
canisters to meet thermal energy storage requirements. 
Detailed heat transfer analyses of the canister are 
necessary to determine temperature histories for subsequent 
use in thermal-stress and material durability calculations. 

The nature of canister heat transfer is very complex. 
Solid-liquid phase change along with many modes of heat 
transfer are encompassed in this time-dependent, three- 
dimensional problem. Although there are several methods 
available for solving classical moving boundary or Stefan 
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problems, a weak numerical solution technique is the only 
feasible approach for the canister problem. This is due to 
the combination of canister geometry and periodic boundary 
conditions which can create multiple, complex-shape phase 
boundaries whose locations are not known a priori. 

Moreover, the enthalpy formulation appears to be the best 
suited weak solution technique to employ on the basis 
accuracy and reliability. 

The TES canister problem has been analyzed by several 
researchers using a variety of different approaches. Many 
of these approaches either 1) over simplify the problem by 
ignoring modes of heat transfer, void effects, and/or free 
convection effects or 2) over complicate the problem by 
rigorously analyzing nearly all facets of canister heat 
transfer in the three-dimensional domain. The need exists 
for canister analyses that provide a balanced approach which 
captures the salient facets of canister heat transfer in a 
step-by-step fashion and analyzes them with a minimum amount 
of required rigor. With the aim of providing timely and 
accurate solutions useful for engineering purposes, this 
approach, described in Chapter II, is adopted for conducting 
canister analyses in the work presented herein. 

In Chapter III, the governing equations for one- 
dimensional semi-infinite PCM, one-dimensional PCM slab 
canister, and two-dimensional (r,z) canister problem 
geometries are developed. Conservation of energy is 
formulated with enthalpy as the dependent variable which can 
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in turn be related to temperature through a set of 
constitutive equations. Void heat transfer is formulated as 
uncoupled conduction and radiation processes. Void shape 
and location are prescribed while void size is determined 
based on conservation of mass. Liquid PCM free convection 
heat transfer effects are modeled through use of a thermal 
conductivity enhancement factor (i.e., the Nusselt number) 
based on existing empirical correlations. 

Chapter IV contains a discussion of the finite- 
difference, simple explicit numerical solution approach 
selected to solve the conservation of energy equation. This 
approach was selected on the basis of simplicity and 
accuracy. Stability requirements and grid size selection 
analyses are also discussed along with a method employed to 
modify the computational domain to account for PCM expansion 
and contraction. 

In Chapter V, numerical solution accuracy is compared 
with available exact solutions and good agreement exists. 
Furthermore, numerical consistency checks confirm that a 
high degree of computational integrity is present in the 
calculations. Initial analyses on one-dimensional canister 
models show that thermal performance is sensitive to the 
type of boundary conditions employed. In addition, the 
effects of void heat transfer and free convection on 
canister performance are shown to be substantial. Two- 
dimensional canister analyses show that the effects of a 
void and free convection are much less pronounced since a 
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large portion heat transfer takes place via conduction in 
canister walls. Thus, the difference in canister 
performance during ground-based tests, in 1-g, and flight 
operation, in micro-g, are predicted to be only moderate. 

In Chapter VI, major conclusions drawn from the one- 
dimensional analyses and two-dimensional analyses are 
listed. In addition, suggested areas for future work are 
discussed . 



CHAPTER II 


INTRODUCTION 

Solidification heat transfer plays an important role in 
many engineering problems. Casting processes, ice accretion 
on vehicles, cryosurgical procedures, structural design in 
permafrost regions, and advanced residential and commercial 
cooling systems are but a few examples. Solid-to-liquid 
phase change materials (PCM's) have also been incorporated 
into the designs of many thermal control and thermal energy 
storage (TES) systems due to their inherent advantages of 
small operating temperature range and efficient, high 
specific energy storage capability. A review of the 
literature yields many references to theoretical and 
experimental work on such systems (see Blumenberg and 
Weingartner (1988), Tanaka et al. (1989), Torab (1989), and 
Sheffield (1981)). PCM TES systems are particularly well 
suited to solar thermal-electric power conversion systems. 

In this application, the ability to adapt the energy supply 
to the energy demand is essential since terrestrial systems 
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must contend with transient cloud cover conditions and 
spacecraft-based systems must adapt to the intermittent 
solar energy supply provided in low earth orbits with 
substantial eclipse periods. 

Perhaps the most notable spacecraft solar power system 
is the one currently under development for the 
NASA/International Space Station Freedom (SSF). SSF 
electrical power will be generated by photovoltaic solar 
arrays initially and later augmented with Solar Dynamic 
Power Modules (SDPM's). The SDPM, shown conceptually in 
Figure 2.1, employs a concentrator to collect and focus 
solar energy into a cylindrical cavity heat receiver where 
it is converted to thermal energy. A fraction of the 
thermal energy is transferred to a circulating working fluid 
to operate the power conversion unit (PCU) (a Brayton cycle 
heat engine) which generates electrical power. The 
remaining thermal energy melts a eutectic composition LiF- 
CaF 2 Phase Change Material (PCM) contained in multiple 
canisters brazed concentrically around working fluid tubes. 
The working fluid tubes run the length of the heat receiver 
cavity which is shown conceptually in Figure 2.2. A single 
PCM containment canister is shown in Figure 2.3. The PCM 
stores and releases thermal energy by undergoing phase 
change at its critical temperature of 1040 K. This permits 
continuous operation of the heat engine during the 
substantial eclipse periods (up to 36 minutes) of Freedom's 
low earth orbit. The design life requirement for the heat 
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LENGTH 299 CM 


Figure 2.2. Heat Receiver. 
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Figure 2.3. PCM Containment Canister. 
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receiver is 30 years. 

A detailed understanding of containment canister heat 
transfer is important to ensure that an efficient heat 
receiver design (that meets all requirements) is developed. 
Three primary technical requirements driving the receiver 
design are: 1) supplying the required thermal power to the 

heat engine working fluid, 2) storing adequate thermal 
energy for use during the eclipse portion of Freedom's 
orbit, and 3) meeting a 30 year design life. The first two 
requirements can be addressed by relatively coarse receiver 
heat transfer analyses. However, the last requirement, a 30 
year receiver design life, is probably driven by canister 
life. Therefore, detailed canister heat transfer analyses 
are required to accurately determine temperature histories 
for subsequent use in thermal-stress and material durability 
calculations. In addition, it is likely that analytical 
models, verified with ground-based experiments, will be 
required to predict on-orbit, flight performance due to the 
limited availability of funds to perform flight experiments. 

To address the need for detailed TES canister analyses, 
a numerical heat transfer model was developed. This thesis 
documents the step-by-step development of this TES canister 
heat transfer model and discusses the numerical results from 
analyses conducted with this model. 
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?.i Attributes of Canist er Heat Transfer 
2.1.1 Thermal Loading 

During nominal TES charge-discharge operation, energy 
is added or removed from the canister outer peripheral 
surface via radiation exchange within the heat receiver 
cavity. The magnitude and sign of this energy exchange 
varies with time and circumferential canister position. 
Energy is removed from the canister inner peripheral surface 
via forced convection cooling by the heat engine working 
fluid, a 39.94 molecular weight helium-xenon gas mixture. 

The temperature of the gas varies with time. Canister 
sidewalls are thermally insulated and can be considered 
adiabatic . 

2.1.2 Role of Conduction Within Canister Walls 

For this TES concept, canister walls are required to 
contain the PCM and to act as effective heat transfer fi 
due to poor PCM thermal conductivity, 0.0382 and 0.0170 
W/cm-K for the solid and liquid phases, respectively. 

Energy is distributed radially, axially, and 
circumferentially by conduction within the canister walls. 
This distribution of heat is required to efficiently heat 
the working fluid and melt the PCM in addition to 
controlling canister wall temperature gradients that give 
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rise to thermal stresses. 

2.1.3 Void Behavior 

During the melting process, the PCM expands 
approximately 20% by volume. As a result, during the 
canister PCM fill process, ullage volume must be left in the 
canister to accommodate melting phase change expansion. 

This ullage volume, or void space, is filled with PCM vapor 
and as a function of time, grows and shrinks during PCM 
freezing and melting, respectively. The void shape and 
location within the canister is determined by a combination 
of surface-tension and buoyancy forces. During ground 
operation in 1-g, it is expected that buoyancy forces would 
dominate and the void would be located in the uppermost 
canister volume. However, it has been shown by post-test 
radiographs that voids associated with PCM solidification 
shrinkage can form on the canister bottom for certain 
cooling conditions, Tong et al . (1987). This situation 

resulted from the combination of high PCM wettability, 
canister geometry, and cooling conditions which permitted 
PCM on the bottom to freeze last. 

The exact shape and location of the void in micro- 
gravity has not been guantified as of this writing. 

However, it is believed that the void shape will be 
essentially spherical (to minimize surface free energy) and 
will be located in the region containing the warmest liguid 
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PCM (liquid LiF-CaF 2 data indicate that surface tension 
decreases with increasing temperature which would result in 
a propulsive force to move a freely suspended void from cool 
liquid to hot liquid). 

2.1.4 Void Heat Transfer 

Across the void, energy is transferred by means of 
conduction, convection, radiation, and 

evaporation-condensation. Scoping calculations have shown 
that PCM vapor convection heat transfer is negligible, 
Whichner et al. (1987). Yet at typical canister operating 
temperatures (950 to 1150 K), conduction, radiation and 
evaporation-condensation heat transfer modes can be 
comparable in magnitude. Kerslake and Ibrahim (1990) showed 
in one-dimensional analyses that void vapor conduction and 
radiation heat transfer are of the same order-of-magnitude 
and are highly dependent on void size. 

2.1.5 PCM Radiant Transmission Characteristics 

There is evidence that suggests significant radiant 
heat transfer through the liquid PCM will likely take place. 
Data show that both solid LiF and CaF 2 PCM components have 
optical "windows" in the 0 to 6 micrometer wavelength range 
where highly polished and monocrystalline specimens exhibit 
*95% transmittance at room temperature. For a black body 
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at 1100 K, *80% of the emissive power occurs at wavelengths 
in this 0 to 6 micrometer window. Thus, a good potential 
exists for substantial radiant interchange between interior 
canister walls. However, the "as-cast" PCM has 
polycrystalline structure and is visually opaque which 
suggests that the shorter wavelength portion of the window 
has been "closed". Hence, radiant transfer through liquid 
PCM regions is likely to be more important than in solid PCM 
regions . 

2.1.6 Convection in the PCM Melt 

Convection in the liquid PCM is driven by buoyancy 
forces, thermal-capillary forces, and by PCM phase change 
expansion/contraction at the solid-liquid interface. Under 
1-g conditions, Whichner, et al. (1987) showed that free 
convective flow in a TES canister is dominant over surface 
tension and advective flows. Also predicted was the 
occurrence of peak wall temperatures located 45 to 90 
degrees around the canister circumference from the location 
of peak heat input at the canister bottom. This occurrence 
was attributed to a vortex shedding mechanism within the 
liquid PCM region which created hot liquid vortices rising 
along the canister outer wall. This prediction was later 
qualitatively confirmed in experiments by Tong et al. (1988) 
where measured peak temperatures occurred at a location 45 
degrees from the canister bottom for a portion of the 
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melt-freeze cycle. In another study, Nusselt numbers (Nu) 
in the 4 to 5 range were predicted for a fully molten PCM 
containment canister during ground tests, Kerslake and 
Ibrahim ( 1990 ) . 

In micro-gravity, thermal-capillary flow is the 
dominant mode of convection. This type of flow arises due 
to surface tension variation along the PCM liquid-void 
interface as a result of temperature gradients. Whichner et 
al. (1987) showed that these flows have an 

order-of-magnitude lower velocity than buoyancy flows in 1-g 
and that the flow field is fairly localized around the void. 
Thus their contribution to overall canister heat in micro- 
gravity is expected to be small. In addition, phase change 
driven flows were predicted to be 7 orders-of-magnitude 
smaller than buoyancy flows in 1-g and thus, need not be 
considered in overall canister heat transfer. 


2.2 Methods For Solving Phase Change Prol 


There are several methods available for solving phase 
change problems which fall into the general classification 
of "moving boundary" or "Stefan" type problems. These 
methods of solution fall into basically four different 
categories: exact, approximate, strong numerical, and weak 

numerical. In short, exact solutions are available for only 
a limited number of inherently one-dimensional problems, as 
in Carslaw and Jaeger (1959), and for problems involving 
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determination of multi-dimensional steady state solid-liquid 
interface geometry, as in Siegel (1982) and Siegel (1985). 
Approximate solutions, i.e. embedding methods, are limited 
to at most two-phase region, one-dimensional problems. 

Strong numerical solutions, i.e. Douglas-Gallie method, 
explicitly solve for the solid-liquid interface position and 
are generally limited to two phase one-dimensional problems 
or with difficulty, single phase, two-dimensional problems. 
This is due to the problem formulation which requires 
simultaneous solution of the heat diffusion equation in the 
PCM solid and liquid regions and the PCM solid-liquid 
interface energy balance equation. The interfacial energy 
balance is formulated using temperature gradient and 
velocity terms which are normal to the PCM solid-liquid 
interface. Hence for multi-dimensional geometries, these 
terms must be evaluated via partial derivatives in the 
coordinate directions. This becomes a difficult task to 
accomplish when the PCM solid-liquid interface position and 
geometry are not known a priori. A second complication 
arises when the second derivatives of temperature in the 
heat diffusion equation must be evaluated near interfaces 
and boundaries where temperature gradients are 
discontinuous. The three-point finite-difference 
approximation to the second derivative relies on all three 
points being in the same medium. However, this is not the 
case in the vicinity of the PCM solid-liquid interface or 
near canister walls. Thus, the three-point scheme must be 
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modified in these areas continuously throughout the computer 
simulation to accurately determine temperatures, Springer 
and Olson (1962). 

Weak numerical solutions, such as augmented specific 
heat or the enthalpy methods, eliminate complications of 
strong numerical techniques since knowledge of the PCM 
solid-liquid interface is not required. In the augmented 
specific heat method, an artificially high PCM specific heat 
value is substituted for PCM regions within an arbitrary 
temperature range, AT, near the melting point. The 
augmented specific heat value is defined by the PCM heat of 
fusion, H , divided by AT. The artificially high sensible 
energy storage (or release) that occurs over the specified 
AT approximates heat of fusion energy storage (or release) 
during PCM phase transformation. However, it is not clear 
how to appropriately choose the value of AT. Selection of a 
small value risks jumping over part or all of the AT 
temperature range and hence, not properly accounting for all 
of the heat of fusion energy. Selection of a large AT value 
is not consistent with the physics of solid-liquid phase 
transformation of a eutectic composition mixture which 
occurs at one discrete temperature. 

In contrast, the enthalpy method uses enthalpy or 
energy content as the dependent variable in the conservation 
of energy equation. Unlike temperature, enthalpy is a 
continuous function across the solid, mushy, and liquid PCM 
regions and thus, can be calculated throughout the entire 
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PCM domain without regard to the location and shape of 
solid-liquid interfaces. Once enthalpy distributions are 
determined, phase front location is contained implicitly in 
the solution. The solid phase exists where the specific 
enthalpy (energy per unit mass), e, is less than 0, liquid 
phase exists where H m <e, and the approximate phase front 
position is located in the mushy zone that exists where 
OsesH,,. Thus, the enthalpy formulation lends itself nicely 
to TES canister type phase-change problems where multi- 
dimensional geometries with periodic boundary conditions can 
produce multiple, complex geometry phase fronts whose 
locations are not known a priori. 

The primary disadvantages of the enthalpy formulation 
are that the PCM solid-liquid interface(s) are not clearly 
defined and that the mushy zone model does not strictly 
apply to the phase change process of a eutectic composition 
mixture. The former introduces some uncertainty in 
temperature gradients in the PCM mushy zone where thermal 
conductivity can only be estimated. The latter introduces 
physics into the problem analysis which do not occur in the 
physical phase change process, i.e. an extended two-phase 
zone is assumed to exist instead of a sharp solid-liquid 
interface. However, both of these disadvantages can be 
minimized to acceptable levels for engineering calculations 
by refining the computational grid on which the calculations 
are performed. 

The strong inherent advantages and benign disadvantages 
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discussed above make the enthalpy method the simplest, most 
logical engineering approach for accurately simulating 
multi-phase, multi-dimensional phase change problems without 
prior knowledge of the geometry of the phase front (s). A 
more complete discussion of the various methods and solution 
technigues can be found in Solomon (1986). 

9.3 Literature Review 

The complex nature of canister heat transfer is 
apparent. It is a formidable task to accurately model all 
facets of the canister phase change heat transfer problem. 
Thus, usually a compromise is made between modeling 
complexity and accuracy of results. Generally, researchers 
have concluded that PCM convection, void effects, and 
three-dimensionality are key features to incorporate into 
TES canister phase change material models. Yet as of this 
writing, analytical results from a three-dimensional model 
that describes TES canister phase change heat transfer with 
PCM convection and void effects have not been published. 
Several computer models of varying sophistication have been 
developed (or are currently under development) to analyze 
this type of PCM Thermal Energy Storage (TES) canister. 

These analyses reported in the literature have modeled many 
aspects of the TES canister heat transfer problem. 

The model described by Solomon (1986) is relatively 
straight forward in that it predicts temperature and phase 
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distributions in the PCM based solely on conduction heat 
transfer. PCM container walls and PCM void formation (due 
to density difference in the PCM solid and liquid phases) 
are not modeled. This model is used to determine the 
feasibility and overall performance of a TES device 
comprised of PCM canisters. 

Tong et al. (1988) modeled transient, three-dimensional 
conduction heat transfer using a finite-element technique 
and solved the problem using the commercially available 
general purpose thermal-structural analyzer program MARC. 

The phase change process was modeled using a modified 
specific heat capacity value over a small temperature range 
above the PCM melting point. However, liquid PCM 
convection and radiation across the PCM vapor void were not 
modeled. Consequently, analytical results generated could 
only be roughly correlated with ground-based test data. 

Using a similar approach, Strumpf and Coombs (1988) used the 
ANSYS general purpose thermal-structural analysis program to 
predict TES canister thermal-stress performance in micro-g. 
The short-fall of using general purpose computer programs is 
the inability to change or add software necessary to explore 
various PCM and void heat transfer modeling techniques. 

Sedgwick et al. (1989) modeled the three-dimensional, 
transient heat transfer of a high length-to-diameter ratio 
annular TES canister containing PCM in a matrix of felt 
metal. The model used an implicit, finite difference 
approach with an iterative solution technique to solve the 
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energy equation. The phase change process was modeled by 
employing an artificially high PCM specific heat value over 
a small temperature range above the PCM melting point to 
simulate the latent heat effects. Use of the felt metal 
more or less uniformly distributes PCM void volume and 
eliminates natural convection effects. Thus, the PCM can be 
analytically treated as a homogenous solid thermal conductor 
with effective material properties dependent on the amounts 
of solid PCM, liquid PCM, and felt metal. 

Viterna (1989) modeled transient, two-dimensional PCM 
heat transfer including conduction and convection in the 
PCM. The continuity, momentum, and energy equations were 
simultaneously solved using a finite element technique with 
a Galerkin formulation (method of weighted residuals). The 
phase change process was analyzed using an enthalpy 
formulation of the PCM conservation of energy equation 
combined with a thermodynamic equation of state. Analytical 
predictions were verified using a variety of published 
results from the literature. 

Wichner et al. (1988) modeled two-dimensional (r,0), 
transient canister heat transfer including conduction, 
convection, radiation, and PCM evaporation-condensation. 

The continuity, momentum, and energy equations were 
simultaneously solved using a simple explicit, finite 
difference technique. The phase change process was modeled 
using the enthalpy method and a prescribed PCM vapor void 
behavior was included for both 1-g and micro-g environments. 
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Work has been continuing to extend the canister model to 
three dimensions and improve the void model. Wilson and 
Flanery (1988) describe the analytical formulation of the 
transient, three-dimensional PCM problem. However, no 
results have been published to date. Although such a model 
offers the potential for a very refined solution, its 
practical utility is diminished by large computer memory and 
execution time requirements as well as extensive computer 
code development/check-out requirements. 

Several conclusions can be drawn from these references. 
First, convection and radiation modes of heat transfer are 
important and must be considered. Secondly, TES canister 
heat transfer is strongly three-dimensional due to 
asymmetric boundary conditions and orientation with respect 
to gravity (ground operation only). Thirdly, as more 
fidelity is built into the canister heat transfer model 
(fluid flow and three-dimensionality), the practical utility 
of the computer code rapidly decreases since computer 
storage requirements and execution times start to challenge 
computer system capabilities. In some cases, insufficient 
computer memory space has been the limiting factor in 
conducting three-dimensional analyses. 

The need exists for a "design-oriented" computer model 
with moderate sophistication to analyze a PCM canister. 

Such a model would have moderate computer memory and run 
time requirements yet would be capable of multi-dimensional 
PCM canister analysis including simplified models of void 
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behavior and liquid PCM convection. This type' of model 
could serve as a canister design tool generating detailed 
temperature distributions for use in structural models and 
for validating less detailed heat receiver models. In 
addition, this model could address key questions about 
canister analyses such as: How should void heat transfer be 

modeled? What effect does the void have on canister heat 
transfer? What are the differences in canister heat 
transfer during ground tests (in the presence of free 
convection) and during flight operation under micro-gravity 
conditions? It seems logical that these questions should 
first be addressed by relatively simplified analyses which 
are likely to yield error-free answers in a timely manner. 
Then, if required, important phenomena identified can be 
modeled in greater detail to refine predictions. 

2.4 Thesis Approach 

In keeping with the “design-oriented" philosophy 
discussed above, the primary thrust of the work herein is to 
develop a PCM canister heat transfer computer code with 
low-to-moderate run time and sufficient accuracy to conduct 
design trade-off or optimization studies and the ability to 
answer the questions posed above. The approach to PCM 
canister code development incorporates an incremental build- 
up of code complexity. This allows the resulting analytical 
predictions to be interpreted without ambiguity by 
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comparison with previously verified solutions. 

Initially, one-dimensional models are analyzed. A 
semi-infinite PCM geometry is first analyzed primarily to 
check the accuracy of numerical methods against a limited 
group a exact solutions. Secondary objectives include 
exploring the effects of applied boundary conditions, void 
heat transfer models, and liquid PCM free convection on the 
solutions to classical Stefan problems. 

A PCM slab canister is next analyzed to evaluate the 
heat transfer performance of an idealized TES canister with 
boundary conditions typical of heat receiver operation. 
Initially, void and free convection models are not 
incorporated into the analyses. Results from these analyses 
are compared with the previously verified results for the 
semi-infinite PCM geometry. Then the numerical model is 
modified to determine the impacts of a vapor void and liquid 
free convection on canister heat transfer performance. 

Finally, a two-dimensional (2D(r,z)) PCM canister model 
is analyzed, first without considering void and free 
convection effects, and then later modified to include these 
effects. Results from these analyses are discussed in a 
comparative manner, highlighting significant differences in 
PCM containment canister temperature and phase distributions 
that arise from the presence of a void and/or free 
convection . 

The transient, multi-dimensional PCM canister heat 
transfer is analyzed using the simple explicit, finite 
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difference numerical method. Conduction, convection, and 
radiation modes of heat transfer are included. The PCM vapor 
void model includes a prescribed void shape and location. 

Void heat transfer occurs via uncoupled vapor thermal 
conduction and internal void surface radiation. To limit 
complexity and computational requirements, liquid PCM flow 
analysis is not performed. Instead, an effective liquid PCM 
thermal conductivity is calculated based on an existing 
Nusselt number correlation. The phase change process is 
numerically analyzed using a solution technique based on an 
"enthalpy" formulation of the conservation of energy equation. 

Canister thermal performance for ground-based (in 1-g) 
and orbital flight (in micro-g) operating modes is predicted. 
Two primary differences in canister PCM behavior are 
anticipated as a consequence of these different operating 
modes: 1) the magnitude and direction of PCM liquid 

velocities and 2) the location and shape of the vapor void. 
During flight operation, it is assumed that only conduction 
heat transfer takes place in the solid and liquid PCM. During 
ground-based operation, it is assumed that conduction heat 
transfer takes place in the solid PCM and that conduction and 
free convection heat transfer take place in the liquid PCM. 

The void shape and location are assumed to be the same for 
ground-based and flight canister operating modes. In both 
operating modes, the void is conservatively located adjacent 
to the canister surface where heat input is applied. 



CHAPTER III 


PROBLEM FORMULATION 


3.1 Problem Statement 


The problem considered in this work is to analytically 
predict the transient temperatures, heat transfer rates, and 
PCM phase distributions in a TES canister comprised of a 
metallic shell containing a eutectic composition LiF-CaF 2 
PCM. The temperatures of the TES canister gaseous cooling 
fluid are also predicted. Conduction heat transfer is 
analyzed in the container walls, solid PCM, and liquid PCM. 
Conduction and radiation (subject to diffuse, gray 
assumptions) heat transfer is analyzed in the void region. 
Void shape and position are specified while void growth or 
shrinkage obeys conservation of mass. Liquid PCM free 
convection is modeled using a modified liquid PCM 
conductivity in a conduction heat transfer analysis. The 
selected problem geometries are a one-dimensional, semi- 
infinite PCM, a one-dimensional PCM slab canister of 
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infinite cross section, and a two-dimensional (r,z) annular 
canister. Constant material thermophysical properties are 
used . 

Phenomena that are not analyzed include PCM vapor 
evaporation-condensation, liquid PCM circulation patterns 
arising from buoyancy or surface tension forces, dynamic 
void shape and position, PCM solid-liquid interface 
kinetics, liquid PCM supercooling, and radiant transmission 

through the PCM. 


•3.2 Governing Eq uations 


3.2.1 PCM Canister Energy Balance 


PCM and canister wall energy redistribution are 
formulated using "the enthalpy method" described by Whichner 
et al. (1988) and Solomon (1986). Based on conservation of 
energy, the governing equation is 

oe ^ = div (kVT) . (3,1) 

dt 

In this equation, e is the specific enthalpy (i.e., given 
Joules per gram), T is the temperature, p is the PCM or 
canister wall density, k is the PCM or canister wall thermal 
conductivity, and t is time. For a special case examined 
with the one-dimensional, semi-infinite geometry, the solid 
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PCM region is translated at a velocity, u, equal to the rate 
of void growth. Therefore, a transport term, u*V(pe), must 
be added to the left hand side of equation (3.1) when 
evaluating the energy balance in the solid PCM. 

When evaluating the discrete form of equation (3.1) in 
the vicinity of the PCM solid-liquid interface or near 
canister wall boundaries, special consideration must be 
given to the conduction heat transfer between adjacent 
finite-difference elements possessing different 
conductivities. This is accomplished by evaluating k in 
equation (3.1) as a "net conductivity", k net , which is 
defined below for the example case of two different material 
slabs placed together in perfect thermal contact: 

k net = V k 2*( s l + s 2 )/( k 1* S 2 +k 2* S 1> * (3,2) 

The net conductivity is based on the individual material 
conductivities, k 1 and k 2 , and conduction path lengths, s 1 
and s 2 . In equation (3.2), materials 1 and 2 could be any 
combination of solid, liquid, or mushy PCM or canister wall 
material . 

As a simplifying assumption, internal PCM radiation 
terms were not included in the solid or liquid PCM energy 
balances. Data indicate that highly-polished and 
monocrystalline specimens of solid LiF and CaF 2 (and 
presumably liquid LiF-CaF 2 ) are semi-transparent to radiant 
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energy with wavelengths less than 6 micrometers. Thus, a 
finite portion of energy is transmitted, absorbed, and re 
emitted within the PCM. 


3.2.2 Constitutive Relationships 

Specific enthalpy is coupled to temperature through 
the following set of constitutive equations: 


’ T . + e / c s 
T = I T m 

' T ,« + ( e - H J/ C L 
, T „ + e / c u 

Here, T is the PCM melting temperature, H B is the PCM heat 
of fusion, and c s , c L , and c w are the specific heat values 
for the solid PCM, liquid PCM, and canister wall material, 
respectively. 


0<e£H Mushy PCM (3.3) 

H <e " Liquid PCM 
_oo<e<<» Canister Walls . 


3.2.3 Mushy Zone Properties 

A so-called "mushy" zone exists when 0<e<H m . This zone 
usually consists of dendritic solid phase surrounded by 
liquid although the exact mushy zone characteristics are 
functions of material properties, temperature gradients, and 
interface kinetics, Flemings (1974) and Grodzka et al. 
(1968). Since extended mushy zones do not exist for 
eutectic mixtures undergoing phase change 
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(i.e., solidification interfaces remain planar or stable to 
within the distance of interlamellar spacing for low 
freezing rates typical of TES PCM's), the mushy zone model 
is only an approximation to the actual phase change process. 
However, this approach greatly simplifies the numerical 
solution technique and the option to shrink the finite- 
difference control volumes to an arbitrarily small size 
(within the limits imposed by computational requirements) is 
available. This, in turn, reduces the mushy zone size and 
hence, reduces the extent of approximation introduced. 

For the purpose of this analysis, the density and 
thermal conductivity of control volumes in the mushy zone 
are treated as linear functions of the liquid PCM volume 
fraction, YF, and mass fraction, XF, such that 


p = ( 1 -YF ) * p $ + YF*p L , 

(3.4) 

k = (l-XF)*k s + XF*k L , 

(3.5) 

where XF and YF are defined as 


XF = e/H, , 

(3.6) 

YF = [ 1 + <p L /p s )*(l/XF - 1) . 

(3.7) 

In these equations, the subscripts S and L denote 
solid and liquid phases, respectively. 

the PCM 
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3.2.4 Void Models 

3.2.4. 1 One-Dimensional Analyses 

For the one-dimensional PCM slab geometry (see Figure 
3.1), the fraction of total canister volume occupied by the 
void, defined as the void volume fraction (WF), varies 
between 0.0 percent, when the PCM is completely liquid, and 
15.44 percent, when the PCM is totally solid. The same 
situation exists for the semi-infinite PCM geometry if an 
arbitrarily large control volume of PCM (or imaginary 
"container") is defined. The PCM growth and shrinkage 
associated with phase transformation is accommodated 
numerically by the combination of variable grid size and a 
variable PCM computational domain. This procedure, known as 
the "combined grid element technique," is described in 
section 4.4. 

Void heat transfer is formulated as conduction, 
radiation, or conduction plus radiation processes. The void 
is assumed to be filled with LiF vapor with negligible 
thermal capacitance and at a pressure equal to the vapor 
pressure of LiF at 1040 K, i.e. 7xl0' 3 torr. These 
assumptions seem reasonable since the vapor pressure of CaF 2 
at 1040 K, as reported by Borucka (1975), is ten orders-of- 
magnitude lower than that of LiF and the void vapor mass is 
very small (10' 8 g) . The void occupies the prescribed 
regions 0 < x £ X v (t) and 6 W £ x s X v (t) for the semi- 
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0 X v x m 

(a) Semi-infinite PCM. 


WALL MUSHY WALL COOLING 

1 VOID LIQUID SOLID 2 FLUID 



(b) PCM slab canister. 


5yy s 0.1 5 cm 
L' =1.30 cm 


Figure 3.1. Schematic One-Dimensional Problem Geometries. 


(a) Semi-infinite PCM 

(b) PCM Slab Canister 
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infinite PCM and slab PCM geometries, respectively. Here 
X v (t) represents the time dependent location of the void-PCM 
interface and 6 y represents the thickness of the PCM 
containment canister wall. 

The time dependent void heat flux, q v , is given by 


q v ( t ) = V 1 *! T(0,t)-T(X v (t),t) ] , (3.aa) 

for the semi-infinite PCM geometry and by 

q v (t ) = R v -1 * [ T(6 u ,t)-T(X v (t) ,t) ] , (3.8b) 

for the PCM slab geometry where R v is the void thermal 
resistance. The void thermal resistance is comprised of two 
components: one associated with heat conduction and one 

associated with radiation. The conduction component of 
thermal resistance is given by X v (t)/k v and [X v (t)-6 g ]/k v for 
the semi-infinite and slab geometries, respectively, where 
the void thermal conductivity, k v , is equal to the thermal 
conductivity of LiF vapor, k LiFVap . Using the kinetic theory 
of gases as done by Wichner et al. (1988), the value of 
k,. BU is 4 . 7x1 O' 4 W/cm-K at 1040 K. 

The radiation component of void thermal resistance, 
assuming gray optical properties and that LiF vapor is a 
non-participating medium, is given in terms of a "radiation 
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conductivity, k rad , by 


x v (t)/k rad = LLZ gp + 1 /„S PCM l]*[T(Q,t)-T(X v m,Ul , (3.9a) 

o*[ T 4 (0,t)-T 4 (X v (t),t)] 


for the semi-infinite PCM geometry and by 


[X v (t)-6 u ]/k rad 




[T(S M ,t)-T(X v (tl 

T*(X v (t),t)] 


m / 


(3.9b) 


for the slab PCM geometry where a is the Stefan-Boltzraann 
constant, € pc(1 is the PCM emittance, and € 0 and € g are the 
emittance values of the surface at x=0 and the containment 
wall at x=6 u , respectively. 

Void heat transfer can be evaluated based on the 
individual conduction and radiation thermal resistance 
components alone or on the basis of an uncoupled, effective 
thermal resistance term incorporating both conduction and 
radiation. Since the void heat transfer components are 
uncoupled, superposition is possible. Using the rule of 
parallel resistances, the effective void thermal resistance 
from conduction and radiation, Ry Eff , is given by 


EFF 



+ . a*(T 4 (0.t>-T 4 (X J (t),tU 

[l/€ 0 +l/€ pc „-l]*[T(0,t)-T(X v (t) 



(3.10a) 


for the semi-infinite PCM geometry and by 
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Rv 


EFF 


‘‘LiFVap 


/[X v (t)-6 W ] 
q * r T 4 ( 6. , , t ) 


[l/€ w +l/€ pcH -l]*[T(« w ,t)-T(X v (t) f t)] 




(3.10b) 


for the slab PCM geometry. Note that the conduction 
component of void thermal resistance is dependent on void 
size and independent of temperature while the converse is 
true for the radiation component. It is also worth noting 
that if the void boundaries of interest consist of PCM only, 
void heat transfer by evaporation/condensation can be 
significant. Scoping calculations by Wichner et al. (1988) 
show that under certain conditions, void heat transfer by 
radiation and vaporization in a LiF PCM are comparable in 
magnitude while heat transfer by conduction is an order-of- 
magnitude smaller. 


3. 2. 4. 2 Two-Dimensional Analyses 

For two-dimensional canister analyses, the WF varies 
between 8 percent, when all PCM is liquid at the melting 
point (T m ), to 22 percent, when all PCM is solid at T B . 

This WF range is the result of receiver fabrication 
requirements and PCM contraction during solidification. 
Unlike the idealized one-dimensional models with no 
additional WF margin, a fraction of the two-dimensional 
canister model volume must consist of PCM vapor void at all 
times during the orbital cycle. The small volume changes 
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associated with cyclic thermal expansion of the PCM and 
containment canister walls are ignored. 

The void geometry selected is a cylindrical annulus 
which easily conforms to a cylindrical finite-difference 
element grid network. The void is placed adjacent to the 
canister outer wall, a location that generates 
conservatively high canister wall temperature distributions 
(see Figure 3.2). Void growth or shrinkage occur uniformly 
across the PCM-void interface defined as r v . As PCM 
liquifies or freezes, r v increases or decreases, 
respectively, about 0.1 cm which changes void volume. An 
attempt was made to accommodate PCM growth-shrinkage in the 
two-dimensional canister analyses by applying a modified 
version of the one-dimensional combined grid element 
technique. However, problems with PCM mass and energy 
balances were encountered. Non-uniform PCM-void interface 
growth-shrinkage approaches were considered, but numerical 
implementation of such approaches are considered beyond the 
scope of the current work. Therefore, as an engineering 
approximation, a constant 15 percent WF is assumed. 

Uncoupled void vapor conduction and radiation are 
considered in two-dimensional canister analyses. Since void 
vapor mass is negligible, the void vapor temperature 
distribution is determined by the steady state heat 
diffusion equation: 


1 1 /r 2 t\ + ( 21 \ 

r dr y dr J dz \dz J 


0 


(3.11) 
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2.261 cm 
1.022 cm 
2.540 cm 
0.152 cm 
0.165 cm 
0.152 cm 


Figure 3.2. Schematic Two-Dimensional Problem Geometry. 
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Preliminary void conduction calculations show that axial 
temperature gradients are small and can be ignored for 
engineering calculations. This eliminates the second term 
in equation (3.11) and the resulting solution has the 
familiar logarithmic form 

T(r) = A*ln ( r ) + B , (3.12) 

with the term A=[T(r 0 )-T(r v ) ]/ln(r 0 /r v ) and the term B=T(r 0 )- 
ln(r 0 )[T(r 0 )-T(r w ) ]/ln(r 0 /r w ) . Equation (3.12) is evaluated 
as a function of time at each axial void grid element to 
determine the void vapor temperature distributions. 

Void radiation heat transfer is calculated based on the 
assumptions that, 1) all void surfaces are diffuse and gray, 
2) PCM surfaces are opaque to all wavelengths of radiation, 
and 3) void vapor is a non-participating medium. With these 
assumptions, the governing equation set for void radiation 
heat transfer, found in Siegel and Howell (1981), is: 


HRS 



NRS 

' 2 F kj o(T k ‘-Tj‘) 

j-1 


(3.13) 


where k indexes from 1 to NRS. In equation (3.13), the 
subscripts k and j are void surface element numbers that 
take on all integer values between 1 and NRS. NRS is the 
total number of radiating surfaces in the void enclosure. 
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The terra A k j is the Kronecker delta equal to 1 for k=j and 
equal to 0 for k*j . Given that each surface element 
temperature, T^, emittance, €^, and element-to-element view 
factor, F kj , are known, the surface element net radiative 
heat loss matrix, Q j /A j , can be determined. These heat loss 
terms are then added to the enerqy balance equation 
(equation 3.1) for the appropriate finite-difference grid 
elements in the canister outer wall, canister side walls, 
and outermost PCM. 

An emittance value of 0.52 is selected for canister 
walls which are fabricated of Haynes alloy 188 (HA 188). 

This value is based on experimental measurements from 
diffuse (grit blasted) HA 188 test coupons for the 
temperature range 1000 K to 1100 K. An emittance value of 
0.6 is selected for PCM surfaces. This value is an estimate 
based on emittance data for similar dielectric materials in 
the temperature range of interest. 

Element-to-element view factors, F k j, were determined 
using existing closed-form view factor solutions and 
considerable view factor algebra (see Rea (1975), Minning 
(1970), Leuenberger and Person (1956), and Sparrow et al. 
(1962)). View factors are recalculated for the various void 
sizes encountered during a simulated melt-freeze cycle. In 
the current work, however, view factors are calculated only 
once for the single void size assumed, i.e. the 15 percent 


WF case. 
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3.2.3 Liquid PCM Free Convection Models 

A rigorous treatment of liquid PCM free convection 
requires simultaneous solution of the three conservation 
equations: namely, conservation of mass, momentum, and 
energy. In light of the presence of solid-liquid phase 
change, void vapor regions, and multi-dimen6ional problem 
geometry, the numerical solution of the TES canister problem 
with PCM liquid circulation becomes extremely complex. 

Since wall temperature distributions (which determine heat 
transfer rates and thermal stresses) are of primary concern 
in TES canister analyses, it seems reasonable to approximate 
the gross behavior of liquid PCM circulation in terms of its 
overall contribution to heat transfer. To this end, a 
substantial simplification in the problem formulation and 
numerical solution is achieved when the conservation of mass 
and momentum equations are eliminated and a modified 
conservation of energy equation is used. 

In the conservation of energy equation, the thermal 
conductivity term, k, is modified based on a simplified 
model of liquid PCM free convection. Free convection models 
are based on existing empirical heat transfer correlations 
in the literature. Enhanced heat transfer due to liquid PCM 
circulation is accounted for by modifying or enhancing the 
value of liquid PCM thermal conductivity, k L , such that 

k LE = Nu*k L 


/ 


(3.14) 
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where k LE is the enhanced PCM liquid conductivity and Nu is 
the Nusselt number. This approach was successfully used by 
Humphries (1974) to predict PCM melt zone height during 
ground testing of a finned thermal capacitor and by other 
researchers (see Szekeley and Chhabra (1970) and Chiesa and 
Guthrie (1974)) to study phase change processes in metals 
and alloy systems. 

There are three concerns to consider with this approach: 
1) exact temperature distributions in the liquid PCM are not 
predicted, 2) empirical Nu correlations for the exact, time- 
varying liquid region geometries and boundary conditions do 
not exist, and 3) existing empirical Nu correlations were 
generated without the presence of phase change. The first 
concern is not critical for these analyses since canister 
wall temperatures are primarily controlled by overall heat 
transfer rates and the solid-liquid interface position and 
not local liquid temperature gradients. The second concern 
also does not appear to introduce major difficulties into 
the analysis. This is based on numerical evaluation of 
several existing correlations shown in Table I which 
indicate that calculated Nu numbers are not extremely 
sensitive to geometry or type of boundary condition. For 
values representative of a canister with fully liquified 
PCM, i.e. Ra=2 . 7* 10 5 and Prandtl number, Pr=2.4, the 
variation in calculated Nu number is only ±13 percent for a 
variety of geometries and boundary conditions. 



Table I. Various Nu Correlations with Ra=2.7xl0 5 and Pr= 
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The third concern has been addressed in experimental 
studies and again, does not appear to present a problem for 
this approach. Kemink and Sparrow (1981) found that 
standard free convection correlations could be accurately 
applied to the problem of PCM melting in open or closed 
containers. In addition, Sparrow et al. (1978) found that 
calculated film coefficients during melting experiments were 
within 12 percent of those calculated for free convection 
experiments with the liquid phase alone. These experiments 
were conducted with the same materials and test apparatus so 
that a direct, quantitative comparison could be made. 

Therefore, the simplest Nu number correlations, for 
horizontal and vertical layers from Ozisik (1985), are 
selected for one- and two-dimensional canister analysis. 
Since the assumption of axisymmetry in two-dimensional 
analyses requires the gravity vector to be parallel with the 
canister axis of symmetry, only the vertical layer 
correlation is used. This restricts the simulated canister 
ground-test orientation to one with the axis of symmetry 
vertical . 

For the semi-infinite PCM geometry, the Nu number 
correlation from Table I for a horizontal layer heated 
isothermally from the bottom is used, Ozisik (1985) . The 
correlation has the form 


n3 


Nu 


C5*Ra' 




(3.15) 
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which is valid for the Prandtl number range 1 < Pr < 20. 
Values for C5 and n3 are given in Table II. Here the 
Rayleigh number, Ra, is defined by 

Ra = q*fl*fT(0.t)-T,1*X , 3 (t) , (3.16) 

a*v 

where g is gravitational acceleration, X (l (t) is the PCM 
liquid zone height equal to the characteristic length, and 
a, (3, v, and T m are the PCM thermal diffusivity, volumetric 
thermal expansion coefficient, kinematic viscosity, and 
melting temperature, respectively. 

The Nu number correlation used for the one-dimensional 
PCM slab and two-dimensional annulus canister geometries, 
also from Ozisik (1985) in Table I, is valid for a vertical 
layer with isothermal or isoflux heating from one side. It 
is given as 

Nu = C5*Ra n3 *(LH/w)" 03 , (3.17) 

with the restrictions of 1 < Pr < 20,000 and vertical layer 
height to width ratio, LH/w, 10 < LH/w < 40. Values of the 
constants are given in Table II. The Ra number in equation 
(3.17) is given as 

Ra = g*/3*[T(X v (t),t)-TJ * (X B (tXz2t v (t)) 3 , 

a*v 


(3.18) 



Table II. Selected Nusselt Number Correlations For Liquid PCM Free Convection 


5 
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for PCM slab canister analyses. Here, the characteristic 
length is the liquid PCM layer thickness which can be 
obtained from the difference between the PCM solid-liquid 
interface and the PCM-void interface, X B (t )-X v (t) . For the 
two-dimensional canister analyses, the Ra number is given as 


Ra = g*/3*[T(r v 



/ 


a*v 


(3.19) 


where the characteristic length, r v -r |> (z,t), is the radial 
liquid PCM layer thickness which is a function of axial 
position, z. The axial dependence of the Ra number is 
removed by substituting integrated average values for the 
PCM-void interface temperature, T(r y ,z,t), and the PCM 
solid-liquid interface position, r m (z,t). 

3.2.6 Canister Cooling Fluid Heat Transfer 

A constant film coefficient, h, is determined based on a 
Nu number correlation discussed by Taylor et al. (1988) 
which is valid for fully developed turbulent flow in 
circular tubes with a low Pr number fluid. The canister 
cooling fluid (or heat engine working fluid), a 39.94 
molecular weight helium-xenon (He/Xe) mixture, has a Pr of 
0.24. The correlation has the form 

= 0 . 022*Re°- 8 *Pr°- 6 , 


Nu 


(3.20) 
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where Re is the Reynolds number. The film coefficient can 
then be evaluated by 

h = Nu*k f /D , (3.21) 

where k f is the cooling fluid conductivity and D is the 
cooling fluid tube inner diameter* 

The cooling fluid mean temperature profile, T f (z,t), is 
determined by a quasi-steady state analysis. This approach 
eliminates extremely small time steps required for transient 
numerical temperature solutions of the cooling fluid which 
has negligible thermal inertia. T f (z,t) is evaluated as a 
function of time such that 

L 

J{U*rr*D*[T w (z, t)-T f (z,t) ] }dz = m*c f * [T f (L, t ) -T f ( 0 , t ) ] , (3.22) 

0 

where T w (z,t) is the canister inner wall-cooling fluid tube 
central temperature, m and c f are the cooling fluid mass 
flow rate and specific heat, respectively, and U is the 
overall heat transfer coefficient given by 

U = [ 1/h + D*ln(l+6 i /D)/(2*k w ) ] -1 . (3.23) 

In this equation, 6 - is the cooling fluid tube plus canister 
inner wall combined thickness and k w is the tube/canister 


thermal conductivity. 
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3.3 Boundary and Initial Conditions 

The boundary and initial conditions for the semi- 
infinite PCM, slab PCM canister, and the two-dimensional 
(r, z ) PCM canister problems are contained in Table III. 

Exact solutions for the Stefan problem are available in 
Solomon (1979) and Solomon (1981) for semi-infinite PCM 
geometries initially at uniform temperature with an imposed 
constant temperature at one face, i.e. problem numbers 1 and 
2 in Table III. A specialized exact solution to the Stefan 
problem with void formation is given in Solomon et al. 

(1986) . 

For the semi-infinite PCM problems 1 and 2 in Table III, 
the Stefan number (St), defined by c*AT/H (| , is selected to 
be 0.10. Here AT is the absolute value of the difference 
between initial PCM temperature and the imposed temperature 
at x=0 . This small Stefan number value is representative of 
phase change processes in TES canisters. For problem 3, the 
value of g is chosen such that the same amount of PCM energy 
change occurs as with the constant temperature boundary 
condition phase change process. 

For the PCM slab geometry, problem 4, values typical of 
a LiF-CaF 2 filled TES canister are selected for the length, 
L’, initial temperature, T s , heat flux input, q(t), film 

coefficient, h, and cooling fluid temperature, T^. Values 
selected also permitted full PCM melting and freezing during 



Table III. Boundary And Initial Conditions 
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the canister simulated charge-discharge cycle. 

Boundary condition values for the two-dimensional 
canister geometry, problem 5, are based on results obtained 
from the heat receiver analysis computer code described by 
Strumpf and Coombs (1988). Figure 3.3 shows the absorbed 
heat input function, q(t), applied to the canister outer 
surface at r=r 0 and the cooling fluid inlet temperature 
function, T f (0,t), for a canister located about 115 cm 
behind the conceptual heat receiver aperture plane (see 
Figure 2.2). During a simulated 91 minute Space Station 
Freedom orbit, canisters in this region of the receiver 
experience maximum heat input and undergo complete PCM 
melting and freezing. Note that q(t) is negative for about 
the first half of the eclipse period when the hottest 
canisters lose heat to the remaining canisters in the 
receiver cavity which are at a cooler average temperature. 
For the second half of the eclipse period, however, 
relatively cold heat engine working fluid from the receiver 
inlet manifold preferentially cools what were the hottest 
canisters at the beginning of eclipse to a temperature below 
that of the remaining canisters in the receiver cavity. 
Therefore, q(t) is positive during this period. It will be 
shown later that the period of canister outer wall heat loss 
for the first half of eclipse significantly affects PCM 
freezing patterns and canister temperature distributions. 

Variations in the inlet cooling fluid temperature, 
T f (0,t), are the result of variations in the average 



2 ujo/m ‘(j)b 
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Sun >|<— Eclipse --->| 



Figure 3.3. Boundary Conditions For Two-Dimensional Analyses. 


T f (0,t), K 
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receiver cavity temperature. It will be shown later that 
interesting local fluid heat transfer rates occur as the 
result of different temperature transients between the 
canister, driven by local absorbed heat fluxes and local PCM 
quality, and the cooling fluid, driven by the average 
receiver temperature and PCM quality. 

Canister side wall boundaries at z=0 and at z=L are 
treated as adiabatic. This assumption is justified by the 
fact that adjacent canisters on a given working fluid tube 
experience nearly the same heating and cooling boundary 
conditions and thus, operate at nearly identical 
temperatures. Furthermore, the canisters are physically 
separated by ceramic paper spacers during tube assembly 
which minimizes any axial heat transfer that could occur due 
to small side wall temperature differences in adjacent 
canisters . 


3.4 Thermoohvsical Properties 

For the purpose of this analysis, constant material 
properties are assumed. These properties are given in Table 
IV. Note that the LiF-CaF 2 PCM and Haynes alloy 188 
(HA 188) containment canister material properties are 
evaluated at 1040 K while the He/Xe working fluid properties 
are evaluated at 900 K. During cyclic canister operation, 
temperatures generally remain within a range ±100 K from the 

Over this limited temperature range, HA 


PCM melting point. 



Table IV. Material Properties 
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39.94 mol wt 
He/Xe 
at 900 K 
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188 thermophysical properties vary by less than 5 
percent while PCM properties generally vary by less than 15 
percent. The only deviation from "flat" property 
temperature dependence occurs for the liquid PCM density 
which decreases by 31 percent over the temperature range 
from 1040 K to 1140 K. However, temperature-dependent 
properties can readily be incorporated into future analyses 
if deemed necessary. 



CHAPTER IV 


NUMERICAL APPROACH 


4.1 Solutioi 


ilqoril 


The simple explicit numerical method is implemented to 
solve equation 3.1. This method is selected primarily 
because of the ease in numerical equation development and 
programming. In addition, Thibault (1985) ranked the simple 
explicit method third best numerical scheme to solve the 
three-dimensional heat diffusion equation in a 
parallelepiped. In this study, nine different numerical 
methods, including explicit, fully implicit, alternating- 
direct ion- implicit (ADI), and Crank-Nicolson methods, are 
compared on the basis of accuracy, ease of programming, 
computation time and computer storage requirements. Ranking 
first and second best are two similar ADI methods which make 
use of the efficient Thomas algorithm to solve the 
tridiagonal system of equations successively in each 
coordinate direction. 
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4.2 Stabil ity Requirements 

Based on the local maximum principle discussed by Solomon 
et al. (1986), the simple explicit scheme is stable as long 
as the time step, At, is chosen such that 

At < Ax 2 / (2 *a) , (4.1) 

for central PCM grid elements in the one-dimensional PCM 
canister analyses. For two-dimensional canister analyses, 
stability is ensured as long as 

At < , Ar 2 / ( 2*ff ul < 4 * 2 > 

*+Ar 2 *k s / ( k s *6 us 2 +k u *6 gs *Az ) 

for canister sidewall grid elements. In equation (4.1), Ax 
is the grid size and a is the PCM thermal diffusivity of the 
solid or liquid phase. This dictates that At values less 
than 0.0375 seconds and 0.0199 seconds must be selected to 
ensure stability for one-dimensional cases with and without 
the presence of free convection, respectively. In equation 
(4.2), Ar and Az are the radial and axial grid spacings and 
ky, 6 US , and a g are the canister side wall element thermal 
conductivity, thickness, and thermal diffusivity. The term 
^ is given by 

41 = Ar [ 1 + 1 

2*r. I ln(r- +1 /r i ) ln(r i /r.. 1 ) 


(4.3) 
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where r- is the radial coordinate of canister sidewall grid 
element i. This dictates that a At value less than 0.0254 
seconds must be selected to ensure stability. 

The exact At value used as input to the two-dimensional 
canister computer program is 0.0234375 or 3/128 seconds. 

This value is essentially the largest number that is less 
than 0.0254 seconds and belongs to the family of fractions 
defined by 

At = (il)/2 (i2) , (4-4) 

where il and i2 are natural numbers. Fractions in this 
family have the unique ability to be converted from decimal 
to hexadecimal format and vice versa without computer round- 
off error. This measure helps to reduce cumulative 
numerical errors in equations containing At, such as energy 
balances, which can become significant after a large number 
of repeated calculations (over 230,000 time steps are used 
for one TES charge-discharge cycle). 

4.3 Grid Selection 

Grid independence tests were performed using the PCM 
slab canister model to determine the appropriate grid 
size for good solution accuracy and resolution. The 
numerical tests were conducted by selecting a fixed time 
step and evaluating several temperature solutions for 
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increasingly smaller grid spacings. Based on these tests, 
~40 grids per cm of PCM were used in one-dimensional 
analyses where computation times are small. For two- 
dimensional analyses, "20 grids per cm of PCM in the radial 
direction (the primary heat transfer direction) and ~5 grids 
per cm of PCM in the axial direction (the secondary heat 
transfer direction) were selected. The larger, two- 
dimensional grid size essentially maintains the solution 
accuracy of the smaller one-dimensional grid but decreases 
solution resolution in order to limit computational time. 

The two-dimensional finite-difference element model is 
shown in Figure 4.1. For this model, the nominal PCM radial 
and axial grid spacings are 0.05115 cm and 0.27940 cm, 
respectively. Note, however, that the radial grid spacing 
is non-uniform. For analyses with the void model, the size 
of two radial grid spacings is adjusted so that adjacent 
void and PCM element boundaries are coincident with the PCM- 
void interface. The location of the PCM-void interface is a 
function of the total PCM mass and the relative percentages 
of solid and liquid PCM that exist at any given time. 

4.4 Combined Grid Element Technique 

Because of the time-varying void size in the one- 
dimensional PCM slab analyses, the PCM computational domain, 
X v (t )£ x£L'- 6 u , must be continually up-dated throughout the 
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Figure 4.1. Two-Dimensional Finite-Difference Element Model. 
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simulation to prevent computations from being made in the 
void region. This preventative measure is taken because 
void elements are considered to be massless which forces 
numerical computations to quickly go unstable and terminate 
computer program execution. Up-dating the PCM computational 
domain is accomplished by calculating the position of X y (t) 
based on conservation of PCM mass and by implementing a 
"combined grid element" technique. Thi6 technique simply 
combines the element that contains the void-PCM interface, 
iv, with the adjacent PCM element, iv+1, to form one larger 
element, iv’, of width Ax v given by 

Ax v = Ax*(l+MF y ) , (4.5) 

where MF V is the mass fraction PCM contained in element iv 
at any given time. Since Ax v > Ax, this grid space 
adjustment does not affect numerical stability. As the void 
front translates during the simulation, the value iv will 
"jump" at discrete instances of time. Once a jump condition 
has been detected, properties of element iv* are updated for 

the future time step based on the average properties of the 
new elements iv and iv+1 from the current time step. 
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4.5 Computer Resource Requir ements 


Computer programs are coded in FORTRAN 77 and executed 
on an AMDAHL 5860 computer using double precision variables. 
Single precision runs were also made to compare run time 
requirements and accuracy. Table V shows the normalized 
computer execution time requirements for running each 
canister analysis program. Single precision runs reduce 
central processor unit (CPU) time requirements by factors of 
1.1 and 1.6 for the one and two-dimensional models, 
respectively. Generally, single and double precision 
calculations are in agreement to within 1 percent for 
temperature predictions and to within 2 to 3 percent for 
melt front predictions. Addition of the void model to the 
two-dimensional canister computer program increases CPU time 
by 25 percent while addition of the free convection model 
has essentially no impact on required CPU time. See 
Appendix A2 for a discussion of the two-dimensional canister 
analysis computer program NUCAM2DV (Numerical Canister 
Jfodel: Two-Dimensional With Void). 



Table V. Normalized Computer Program Execution Times In CPU Seconds/Simulation Seconds 
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CHAPTER V 

RESULTS AND DISCUSSION 
5.1 Numerical Solution Accuracy 

Numerical solutions were obtained to the one- 
dimensional, temperature controlled freezing and melting 
problems, problems 1 and 2 in Table III, respectively, in 
which exact solutions also existed. Exact and numerical 
temperature and phase front solutions were then compared to 
assess the accuracy of the numerical computations. That 
comparison is shown in Table VI. Numerical melt front 
solutions are within 0.6 percent of the exact solutions 
without voids and within 1.8 percent of the exact solution 
with void. Numerical temperature solutions are within 0.5 K 
of the exact solutions for all problems. These results 
indicate that the numerical scheme is accurate and properly 
implemented . 

To assure proper implementation of the numerical 
equations, the two-dimensional canister computer program 
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Table VI. Numerical Versus Exact Solutions For Problems 1 And 2 
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(without a void model) was exercised independently in the r 
and z coordinate directions and the results were compared to 
the exact solution of problem 2 in Table III. As with the 
one-dimensional analyses, numerical and exact solutions 
agreed to within 0.6 percent. In addition, the computer 
model global energy balance was checked to assure that 
boundary conditions were properly implemented. An energy 
balance was maintained to within 0.003 percent. 

In the absence of applicable exact or analytical 
solutions, previous numerical solutions, and experimental 
data, numerical consistency checks were performed to assess 
the validity of numerical solutions from two-dimensional 
analyses with the void model . A numerical check of canister 
model energy balance and void surface element view factor 
summation was carried out for each computer run. For all 
cases, an energy balance was maintained within 0.0015 
percent and all surface element view factors summed to 1.0 
within machine accuracy. 

5.2 One-Dimensio nal Analyses 


5.2.1 Semi-infinite PCM 


5. 2. 1.1 Effects of the Void 


The temperature controlled freezing problem (problem 1 
in Table III) was solved for three cases: 1) without a 
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void, 2) with a void considering only vapor conduction heat 
transfer, and 3) with a void considering only radiation heat 
transfer. The effect of a void on the PCM freezing process 
is evident in Figure 5.1 which illustrates freeze front 
position, X m , versus time for the three cases outlined 
above. Here, the presence of the void reduces the amount of 
PCM frozen by factors of "4 or "5 assuming void heat 
transfer via conduction only or radiation only, 
respectively. For small St number freezing processes, the 
amount of energy removal is essentially proportional to the 
amount of PCM frozen. Therefore, in this problem where void 
size is small (i.e., 15 percent of X m ), the magnitude of 
void heat transfer via conduction and radiation are 
comparable . 

5.2. 1.2 Effects of Boundary Conditions 

The effect of boundary conditions on the freezing 
processes in problems 1 and 3 of Table III is examined next. 
For these problems, comparison of results between constant 
temperature and constant heat flux boundary conditions is 
made on the basis of equal energy removal. This is 
accomplished by first integrating the boundary heat flux at 
x=0 over the 50 minute period for the constant temperature 
case to determine the total energy removed from the PCM. 

This total energy is then divided by 50 minutes to determine 
the required boundary condition at x=0 for the constant 
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cooling heat flux case. This basis of comparison is 
consistent with a thermal energy storage application in 
which a finite amount of energy must be stored and released. 

Plots of boundary temperature, T(0,t), and freeze front 
location, X B , versus time are contained in Figures 5.2 and 
5.3, respectively. Boundary temperature decreases linearly 
versus time with constant wall heat flux (see Figure 5.2). 
However, with radiation void heat transfer, variation in 
T(0,t) is small. This confirms the anticipated 
insensitivity of radiation heat transfer to void size. 

Freeze front advancement differs substantially with boundary 
condition assuming conduction void heat transfer (see Figure 
5.3). A constant temperature boundary condition generates 
PCM freezing “ time 1/2 and a constant flux boundary condition 
generates PCM freezing ~ time. This freezing behavior is 
characteristic of one-phase Stefan problems (see Yao and 
Prusa (1989)). 

5. 2. 1.3 Effects of Free Convection 

The effect of free convection on a melting process with 
St=0.10 (see problem 2 in Table III) is illustrated in 
Figure 5.4. Figure 5.4 contains plots of melt front 
position, X B , versus time, PCM temperature at time 25 
minutes, T(x, 25 min), versus position, and boundary heat 
flux, q(0,t), versus time for cases with and without free 
convection present. Although constant temperature is 
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Figure 5.3. Melt Front Locations For Equal Energy Removal Processes. 
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Figure 5.4. Problem 2 With St=0.10. 

(a) Melt Front Position 

(b) PCM Temperature Distributions 
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(c) Boundary heat flux. 


Figure 5.4. Problem 2 With St=0.10. 
(c) Boundary Heat Flux 
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maintained at x=0, the progression of X„ is nearly linear 
with the presence of free convection. This is 
characteristic of a conduction controlled, constant heat 
addition melting process (see Figure 5.4a). This same 
behavior was observed in experiments by Hale and Viskanta 
(1980). At 25 minutes, temperature gradients in the liquid 
PCM are reduced by a factor of ~2 while at the same time 
q( 0 , 25 min) is increased by a factor of "3 (see Figures 5.4b 

and 5.4c) . 

In Figure 5.4c, a local peak in q(0,t) exists at 5 
minutes with free convection present (see the top curve). 
This is in sharp contrast to the monotonically decreasing 
behavior of q(0,t) predicted when accounting for only liquid 
PCM conduction heat transfer (see bottom curve). This heat 
transfer over-shoot phenomenon was also observed 
experimentally by Sparrow et al. (1978) and predicted 

numerically by Sparrow et al. (1977). 

The over-shoot occurs when the magnitudes of convection 
and conduction in the liquid PCM are equal. The phenomenon 
can be explained physically in the following way. During 
the early stages of melting, liquid PCM heat transfer is 
controlled by conduction and decreases rapidly as the fluid 
layer thickness, or conduction path, increases. Fluid 
velocities are small and liquid temperature profiles are 
nearly linear. As melting proceeds, fluid velocities 
increase and boundary layers start to form. Moderate 
temperature gradients exist in both the bulk fluid and in 
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the developing boundary layers near the heated and cooled 
surfaces. During this stage of melting, the rapid fall-off 
in conduction heat transfer is compensated for by the 
increasing convective heat transport. Thus, overall heat 
transfer rates start to increase. Fluid conduction and 
convection heat transfer components eventually become equal 
in magnitude at which time the liquid PCM heat transfer rate 
is locally maximized. 

As melting proceeds still further, heat transfer rates 
fall-off again, but at a much slower rate consistent which 
the gradually increasing flow resistance associated with the 
growing liquid region size. Boundary layers become fully 
developed and fluid motion and temperature gradients are 
confined to narrow layers adjacent to the heated and cooled 
surfaces. Temperature gradients through the bulk liquid are 
essentially zero as the heat transport through the liquid 
layer is due solely to liquid recirculation. 

A similar situation arises in constant flux melting 
experiments. In these tests, the measured heat source 
temperature exhibits a maxima (i.e., the film coefficient 
exhibits a local minima) near the transition from conduction 
dominated to convection dominated liquid heat transfer, 
Goldstein and Ramsey (1979). 

It is interesting to note for the case with free 
convection, that although the liquid PCM region continues to 
grow in size, q(0,t) reaches a near steady state value by 50 
minutes (see Figure 5.4c). This suggests that the free 
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convection film coefficient and liquid PCM region thermal 
resistance are independent of X m . A check of the free 
convection correlation in Table II reveals that for Ra>10 , 
the Nu number essentially increases linearly with X t . This 
is to be expected since in this Ra number range, Nu Ra 
X. 0 - 9 . Thus, the film coefficient exhibits a very weak 
dependence on the melt zone size during the later stages of 
melting when large Ra numbers exist. Experiments by Bathelt 
et al. (1979) also confirm this behavior. 

5.2. 1.4 Observations 

Although semi-infinite PCM analyses were performed 
primarily to check the accuracy of numerical predictions, a 
few observations regarding the general performance of TES 
canisters can be made. First, dramatically different 
problem solutions are obtained depending on whether a void 
is present or not. Consequently, any reasonable analysis of 
a PCM with appreciable volume change (i.e., >5 percent) must 
include a model of void behavior. Furthermore, the void 
model must properly take into account the primary modes of 
heat transfer for the given void size and void vapor 
thermophysical properties. 

This point is illustrated by considering two separate 
LiF vapor void problems that lead to two separate 
conclusions: 1) the freezing process shown in Figure 5.1 

where void size is small (i.e., < 0.1 cm) and 
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the temperature difference across the void is moderate 
(i.e., *40 K) and 2) scoping calculations by Whichner et al. 
(1988) where a large void size is considered (i.e., 1 cm) 
and a large void temperature difference exists 
(i.e., "100 K) . In the first problem, it is concluded that 
void heat transfer by radiation and conduction are nearly 
equal since nearly equal quantities of PCM are frozen in 
each case. However, in the second problem, radiation heat 
transfer is predicted to be *35 times larger than conduction 
heat transfer. Therefore, void heat transfer modeling must 
be consistent with the void geometry analyzed and the 
anticipated temperature conditions. A further discussion of 
void heat transfer is given in the next section. 

A second observation can be made concerning the wide 
range of solidification rates, boundary temperatures, and 
boundary heat fluxes predicted for equal energy removal 
processes depending on the type of boundary condition 
assumed (see Figures 5.2 and 5.3). This suggests that the 
type of boundary conditions that a TES canister experiences 
will influence, to some extent, items such as solid PCM 
crystalline structure and void distribution (which are both 
functions of the solidification rate). 

PCM containment canister material durability will also 
be influenced by the type of boundary condition since 
results indicate that a constant temperature heat sink (or 
heat source), such as a nearly isothermal heat pipe, 
maintains the boundary temperature closer to T B than does a 
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constant flux boundary condition. The practical implication 
of this result is that lower PCM containment wall 
temperatures will occur with heat pipe receiver concepts 
versus direct absorption receiver concepts. Keeping 
canister wall temperatures close to T B is important for 
achieving long component design lives since increasing 
temperature enhances the PCM corrosion rate and reduces 
canister material strength. However, the benefit of lower 
canister wall temperatures must be weighed against the added 
complexity and mass of a heat pipe receiver concept. 

A third observation can be made regarding the influence 
of free convection. As shown in Figure 5.4, the heat 
transfer differences from liquid PCM circulation (in terms 
of melting rates, temperature gradients, and boundary heat 
fluxes) are substantial. Therefore, PCM analyses must 
include a free convection model to enable correlation with 
ground-based experiments. This same conclusion has been 
reached by researchers referenced in the preceding section 
and by others (see Bathelt et al. (1979), Deal and Solomon 
(1981), and Humphries (1974)) after completing phase change 
experiments which focus on the effects of free convection in 
the PCM melt. Furthermore, it is essential to be able to 
accurately predict ground-based performance of flight-design 
TES units since full scale flight tests may not be practical 
on the basis of cost. If this is the case, extrapolation of 
TES unit flight performance (in micro-g) can be calculated 
with a satisfactory confidence level. A further discussion 
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of this point is given in the following sections. 

5.2.2 PCM Slab Canister 

5.2.2. 1 Void Thermal Resistance 

Void thermal resistance is plotted as a function of time 
in Figure 5.5 for a representative size slab PCM containment 
canister undergoing a 30 minute, constant heat input melting 
period and a 20 minute, zero heat input freezing period. 
Conduction and radiation thermal resistances are of the same 
order-of-magnitude and both vary significantly with time due 
to variations in void size and canister wall temperature, 
respectively. Therefore, for TES canisters of the type 
analyzed herein, void heat transfer is most accurately 
modeled as a simultaneous conduction-radiation process. The 
resultant thermal resistance from uncoupled conduction and 
radiation heat transfer modes is also shown in Figure 5.5. 

5. 2. 2. 2 Wall 1 Temperatures 

The canister wall 1 temperatures, T(0,t), are shown 
versus time in Figure 5.6 for different void heat transfer 
assumptions. Note that canister heat input is applied at 
wall 1 which is adjacent to the void while wall 2 is 
convectively cooled by the heat engine working fluid (see 
Figure 3.1). T(0,t) predictions widely vary depending on 
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Figure 5.5. Void Thermal Resistance. 
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Figure 5.6. Wall 1 Temperatures With Various Modes Of Void Heat Transfer. 
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the assumed type of void heat transfer. With combined mode 
void heat transfer, T(0,t) is nearly isothermal at 1290 K 
during melting and at 1050 K during freezing. Ignoring the 
conduction component of void heat transfer increases 
canister wall 1 temperature predictions by 50 to 150 K over 
those with combined mode void heat transfer. The magnitude 
of this temperature difference is probably not acceptable 
from a canister design view point. Therefore, ignoring void 
vapor conduction is not a good assumption in this case since 
it would lead to an overly conservative canister design. 

Ignoring void radiation results in wall 1 temperatures 
that exceed the melting range of Haynes alloy 188 (1575 to 
1630 K) . In all cases, wall 1 temperature predictions are 
too high for long term operation of containment canisters 
constructed with superalloys. This illustrates the need for 
heat transfer enhancement fins between heat addition and 
heat removal surfaces when dealing with low conductivity 
PCM ' s to maintain maximum wall temperatures below 1150 K. 

5. 2. 2. 3 Effects of Void Distribution and Consequences 

for One-Dimensional Analyses 

Since the void in actual PCM containment canisters is 
not evenly distributed around the circumference, Strumpf and 
Coombs (1988), the behavior of localized canister radial 
segments can roughly be approximated by the behavior of one- 
dimensional PCM models with or without a void. Figure 5.7 
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illustrates wall 1 temperature predictions versus time for 
cases with and without a void. The void increases wall 1 
temperatures between 50 to 200 K throughout the TES charge 
period with constant heat input. This introduces the 
potential that for PCM canisters with asymmetric heat input, 
the maximum wall temperatures will not occur at the point of 
maximum heat input but instead will occur in the region of 
the void (i.e., the location of largest thermal resistance). 
Hence, the position of the void must be quantified to 
accurately characterize wall temperatures of canisters with 
high length-to-diameter ratio, i.e. canisters with small 
side wall end effects. However, accurate prediction of the 
void location in micro-g requires very complex calculations 
which are difficult to verify. 

To avoid the difficulties in predicting void behavior, a 
straight-forward approach could be adopted in which the void 
is placed adjacent to the heat input surface to yield 
conservative temperature predictions. It will be shown in 
the two-dimensional canister analysis sections that for a 
low length-to-diameter ratio canister, i.e. 1/d ** 0.5, wall 
temperature sensitivity to void location is greatly reduced 
because of the large heat transfer contribution of canister 
sidewalls . 

It should be noted that one-dimensional canister 
analysis accentuates wall temperature increases introduced 
by a void since all canister absorbed energy must be 
transferred across the void. In an actual canister, 
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absorbed energy in the outer wall has multiple heat transfer 
paths in which to diffuse. Thus, wall temperature increases 
due to the presence of a void would be much less pronounced 
than indicated from one-dimensional predictions. 

The sensitivity of wall temperatures to void location 
and to the nature of void heat transfer depends on the 
extent to which PCM heat transfer is required for energy 
redistribution within the canister. In section 5.3, two- 
dimensional analytical results will show that during the 
cycle heating period, when highest canister wall 
temperatures exist, roughly 30 to 70 percent of canister 
total radial heat transfer occurs within the side walls. 
Thus, it will be shown that the sensitivity of canister wall 
temperatures to PCM-void distribution is greatly reduced 
over the one-dimensional case. 

It is conceivable that the magnitude of this reduction 
may render void heat transfer secondary in importance to 
canister and PCM heat conduction/convection. In this case, 
wall temperature predictions would be essentially 
independent of the method used to model void heat transfer. 
Results from steady state PCM canister heat transfer 
analyses discussed by Tong et al. (1988) support this 
assertion. In this study, the maximum canister wall 
temperature increased by only 29 K with the addition of a 
circumferential void at the canister outer diameter. A more 
detailed discussion of canister wall temperature 
sensitivity to a void will be given in section 5.3. 
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5. 2. 2. 4 Effects of Free Convection 

Figure 5.8 illustrates the impact of free convection 
on PCM melt and void front positions, canister wall 
temperatures and PCM temperature distributions. The 
PCM occupies the region 0.15 < x < 1.15 cm between 
the two canister walls. The void occupies the region 
0.15 < x < 0.30 cm as it sequentially grows and shrinks (due 
to density differences in the solid and liquid PCM) during 
PCM freezing and melting, respectively. Little difference 
exists between melt/void front locations with the addition 
of free convection (see Figure 5.8a). However, the presence 
of free convection significantly lowers canister wall 1 
temperatures and PCM temperature gradients during the TES 
charge period in addition to melting slightly more PCM (see 
Figures 5.8b and 5.8c). 

Figure 5.9 shows the liquid PCM Nu number as a function 
of melt front position X m with and without a void. Without 
a void, the critical Ra number is exceeded at X m — 0.50 cm and 
the Nu number increases linearly with X^ until complete PCM 
melting has occurred. At this point, Nu— 3.4 which is about 

35 percent lower than the Nu number for the semi-infinite 
PCM geometry with the same X m . This seems reasonable since 
for the semi-infinite PCM geometry, a horizontal liquid 
layer with heat input from the bottom (with respect to 
gravity) is assumed. This orientation creates greater 
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Figure 5.8. PCM Slab Canister Results. 
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convective instability than expected for the PCM slab 
canister which is assumed to have a vertical liquid layer 
heated in a direction normal to the gravity vector. During 
freezing, the Nu number falls off rapidly with decreasing X B 
and becomes linear with a slope about 25 percent lower than 
during melting. This is attributed to smaller Ra numbers 
due to reduced temperature gradients. This suggests that h 
is essentially independent of X m but does depend on whether 

the PCM is melting or freezing. 

Note that in Figure 5.9, there is a marked difference in 
the curves for the cases with and without a void. During 
PCM melting, Nu numbers for the case with a void appear to 
be lower than those occurring without a void for a given 
melt front position, X m . This result has no physical 
significance but instead is the consequence of how X B is 
measured in the case with a void: namely, the value of X m 

is necessarily increased by the size of the void at any 
given time, i.e. by 0.15 cm at the start of melting which 
vanishes to 0.0 cm at the conclusion of melting (see Figure 
3.1b). 

During PCM freezing, the larger Nu numbers in Figure 5.9 
for the case with a void (compared to the case without a 
void) can be ascribed to physical differences between the 
two problems in addition to the convention adopted for 
measuring X m . With a void present, slightly less PCM is 
liquified during the charging period. Instead, this energy 
manifests itself sensibly in the from of substantially 
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increased wall 1 temperatures. When the PCM is discharging, 
however, wall 1 rapidly cools and effectively acts as a heat 
source to the liquid PCM. This sustains liquid PCM 
temperature gradients for a longer period of time and hence, 
Nu numbers greater than 1.0 persist for smaller values of X B 
during freezing with a void as compared to without a void. 

5. 2. 2. 5 Ground-Based Testing of Flight Design Hardware 

Several comments can be made about the effects of liquid 
PCM free convection and its implications for ground-based 
testing of the conceptual heat receiver or PCM containment 
canisters designed for operation in low earth orbit. 

Analyses have shown that liquid PCM convective flows arising 
from buoyancy and surface tension forces are small in a 
micro-gravity environment, Whichner et al. (1988). Thus, 
liquid PCM heat transfer during on-orbit operation will take 
place primarily via thermal conduction. However, for 
ground-based tests, the effects of free convection (based on 
one-dimensional analyses) are lower canister wall 
temperatures and increased PCM melting rate during heat 
input periods. These effects are enhanced for the canister 
orientation in which the direction of outer wall heat input 
is from the bottom (with respect to gravity) as opposed to 
normal to the gravity vector. Therefore, free convection 
effects should lead to improvement in overall receiver PCM 
utilization, greater receiver cavity isothermal lity , and 
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lower receiver heat losses. During thermal discharge 
periods, free convection effects are small and should not 
significantly affect receiver thermal performance. 

These results suggest that canister ground tests should 
be conducted with an orientation that permits outer wall 
heat input from the top (with respect to gravity) or normal 
to the gravity vector to minimize free convection effects. 
Furthermore, outer wall heating in a direction normal to 
gravity places the canister axis of symmetry parallel to the 
gravity vector. In this orientation, free convective 
effects will tend to be more uniform around the canister 
circumference when compared to the canister orientation with 
outer wall heating directed parallel to the gravity vector. 

It is interesting to note that an analogous situation 
arises in adiabatic, two-phase (liquid-vapor) flow in 
circular tubes. Researchers have found that vapor bubble 
shapes and distributions in 1-g, vertical tube flow tests 
closely match those encountered during low-g flow tests 
while horizontal 1-g flow tests generate substantially 
different vapor bubble characteristics, Siegel (1967). This 
observation introduces the possibility that for certain 
cases, vertical orientation testing in 1-g offers an 
adequate test simulation of anticipated micro-g operation. 

The argument for canister ground test orientation can 
also be extended to heat receiver ground testing. The 
preferred heat receiver ground test orientation should be 
with the axis of the receiver vertical (see Figure 2.2). 
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This orientation permits canisters on each working fluid 
tube to receive outer wall heating directed normal to 
gravity and thereby experience uniform free convection 
effects. A test conducted with the receiver axis horizontal 
would cause canisters on top receiver tubes to be heated 
diametrically opposed to gravity (i.e., generating maximum 
convective activity) while canisters on bottom tubes would 
be heated in the direction with the gravity vector (i.e., 
generating minimum convective activity) . This situation 
would skew tube-to-tube canister performance and introduce 
additional receiver cavity circumferential temperature 
variations not expected during on-orbit operation. A more 
detailed discussion of liquid PCM free convection effects 
based on two-dimensional canister analyses will be given in 
section 5.3. 


5.3 Two-Dimensional Analyses 
5.3.1 Canister Without Void or Free Convection Models 
5. 3. 1.1 PCM Phase Distributions 

PCM phase distributions are shown in Figure 4.10 at 
several times (24.28, 54.63, 66.77, and 91.05 minutes) for a 
91 minute cycle in which the PCM is being charged for the 
first ~55 minutes and discharged for the remaining “36 
minutes. The heat transfer benefits of the canister side 
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Figure 5.10. PCM Phase Maps. 

(a) Time=24.28 min., MFL=0.30 

(b) Time=54.63 rain., MFL=0.99 
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Figure 5.10. PCM Phase Maps. 

(c) Time=66.77 min., MFL=0.63 

(d) Time=91.05 min., MFL=0.00 
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walls are evident at 24.28 minutes into the TES charge 
period. The liquid and mushy PCM regions extend radially 
inward adjacent to canister sidewalls to a greater extent 
than the bulk PCM. This indicates the manner in which 
sidewall heat transfer enhances PCM melting in both radial 
and axial directions without the need for large liquid PCM 
temperature gradients. Nearly complete PCM melting occurs 
by 54.63 minutes at which time only a small mushy region 
exists at the canister inner radius. 

At 66.77 minutes, mushy PCM and solid PCM regions 
completely surround a liquid PCM core region. This phase 
distribution is the result of heat removal at both inner and 
outer radial canister surfaces during the first half of the 
TES discharge period (see Figure 3.3). As freezing 
continues, the solid region growths inward from all sides 
until the liquid core region is consumed at about 6 minutes 
prior to the end of the discharge period. Thus, at 91.05 
minutes, only solid PCM exists. In short, the PCM 
solidification process obeys 2 simple rules: 1) solid PCM 
forms on cooled surfaces and 2) the last liquid to solidify 
is situated furthest from cooling surfaces. 

The simple observations of the PCM freezing pattern 
mentioned above have important implications for evaluating 
void behavior. Knowing where solid PCM formations occur 
narrows down the possible locations that voids can occupy. 
The above freezing behavior suggests that had PCM density 
differences been accounted for, the resultant void volume 
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would most likely end up as a central core region surrounded 
by solid PCM. This prediction is based on the assumption 
that liquified PCM can creep into canister corners and 
completely cover available internal canister wall surface 
area. In this sense, the void in completely liquified PCM 
would be centrally located within the canister volume prior 
to the start of PCM freezing. 

Langbein et al. (1990) found from micro-g experiments 
and Concus and Finn (1990) proved mathematically that a 
liquid will creep into container corners if the sum of the 
liquid contact angle plus the corner half angle is less than 
90 degrees. This situation does in fact exist for liquid 
LiF-CaF 2 at temperatures below *1100 K. Furthermore, 
ground-based observations of canister PCM distributions by 
Strumpf and Coombs (1989) show a centrally located void 
position after repeated freeze-thaw cycles in an air furnace 
where cooling takes place on all canister surfaces. In 
other ground-based experiments by Blumenberg and Weingartner 
(1988), LiF-filled coaxial cylinders were found to have 
voids located at or near non-cooled surfaces. In 
experiments by Sparrow et al. (1978), guard heaters were 
incorporated into the apparatus specifically to control void 
formation during PCM solidification in preparation for 
subsequent melting experiments. As intended, voids formed 
near the guard heaters and therefore, essentially eliminated 
solid PCM porosity and potential over-stress conditions in 
the PCM containment vessel from undesirable void 
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distributions. In the present analyses, cooling occurs 
along all canister walls which confirms the above assertion 
that void formation will likely result in the canister 
central volume. 

5.3. 1.2 Temperature Distributions 

Canister temperature contour maps that correspond to the 
PCM phase maps are shown in Figure 5.11. The isotherms just 
above and below 1040 K reveal the approximate position of 
the PCM solid-liquid interface. Noticeable isotherm 
compression occurs in the vicinity of the melt front as 
evidence of the relatively high heat transfer rates needed 
to support PCM melting (see Figure 5.11a) or PCM freezing 
(see Figure 5.11c). At times when only liquid PCM or only 
solid PCM exists, isotherms are spaced in a relatively 
uniform fashion (see Figures 5.11b and 5. lid, respectively). 
Noticeable bending in the isotherms occurs near canister 
walls. This illustrates the effect of canister wall heat 
transfer enhancement that effectively behaves as a heat sink 
for the outer wall and as a heat source for the inner wall. 


5.3. 1.3 Temperature and Heat Transfer Variations 

Figure 5.12 illustrates the variation in maximum 
canister wall temperature and heat transfer to the cooling 
fluid as a function of cycle time. The maximum canister 
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(b) Heat transfer to cooling fluid. 


Figure 5.12. Cyclic Variations In Maximum Canister Wall 
Temperature And Heat Transfer To The Cooling Fluid. 

(a) Maximum Canister Wall Temperature 

(b) Heat Transfer To Cooling Fluid 
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wall temperature occurs at the axial midpoint of the outer 
wall. This temperature is a strong function of the phase 
change process (see Figure 5.12a). As long as liquid and 
solid PCM coexist (in this case from 9-86 minutes), absorbed 
canister energy manifests itself as latent heat and hence, 
temperatures will not strongly deviate from the PCM melting 
temperature, T B =1040 K. Finite temperature deviations from 
1040 K are required to transfer heat to and from the solid- 
liquid interface. However, these deviations are held to 
acceptable levels by choosing a reasonably high canister 
conductance, i.e. by limiting canister size and selecting 
adequate wall thicknesses. Once a single PCM phase exists, 
large temperature transients result as a consequence of 
sensible energy change. 

Variation in cooling fluid heat transfer can generally 
be ascribed to variation in cooling fluid inlet temperature, 
T f ( 0 , t ) (see Figure 5.12b). Heat transfer to the fluid is 
proportional to the temperature difference between the 
cooling fluid tube wall (i.e., the canister inner wall) and 
the fluid. Since for most of the cycle time, two-phase PCM 
exists, tube wall temperatures remain fairly constant near 
T m . Thus during this period, the temporal change in cooling 
fluid heat transfer is inversely proportional to the 
temporal change in T f (0,t): hence, an increase in T f (0,t) 

leads to a corresponding decrease in cooling fluid heat 
transfer. The exception to this behavior occurs at the 
beginning and at the end of the 91 minute cycle when only 
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solid PCM exists. During these periods, tube wall 
temperature transients are larger than transients in 
T f (0,t). Thus at the beginning of the cycle, for example, 
this behavior results in increasing cooling fluid heat 
transfer with increasing T f (0,t). 

5. 3. 1.4 Side Wall Heat Transfer Fractions 

Since the LiF-CaF 2 PCM is a poor thermal conductor, 
highly conducting canister walls are necessary to distribute 
energy absorbed at the canister outer surface to the PCM and 
to the cooling fluid (heat engine working fluid) without 
excessive temperature gradients. One measure of the 
effectiveness in which the canister walls redistribute 
absorbed energy is the fraction of total canister radial 
heat transfer which occurs via the canister side walls. 

This "side wall fraction" is plotted versus time in Figure 
5.13 for three radial locations: r^, r 0 ‘, and (r^ + r 0 )/2. 

Side wall fractions generally run between 40 and 60 percent 
during the heat input period. Three distinctive dips in the 
curves are evident at times 12 minutes, 32 minutes, and 53 
minutes for locations r 0 ~, (r i + r 0 )/2, r^, respectively. 

These dips are associated with the passage of the PCM melt 
front at which time radial PCM heat flow increases (thereby 
decreasing the side wall fraction) to support the melt front 
advancement . 

During the heat removal period, side wall fractions are 
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widely varying. The side wall fraction at r 0 ~ closely 
follows the variation in boundary heat flux, q(t). 

Initially, q(t) is negative which freezes a thin, isothermal 
PCM layer on the canister outer wall while a small, but 
finite, side wall temperature gradient exists. Thus, a side 
wall fraction in the 100 percent range is achieved. At ~62 
minutes, the thin PCM layer has completely solidified 
creating relatively large PCM temperature gradients. At the 
same time, the negative outer wall absorbed heat flux, q(t), 
combined with inner wall convective cooling effectively 
flattens canister side wall radial temperature gradients 
near the outer radius. Thus, the side wall fraction 
approaches zero during this period. At 73 minutes, q(t) 
becomes positive and quickly re-establishes side wall 
temperature gradients forcing the side wall fraction up into 
the 70 to 100 percent range. 

The opposite behavior occurs at the canister mean 
radius, (r i + r 0 )/2, during the heat removal period. As 
shown by Figure 5.10c, a small liquid PCM zone surrounded by 
solid and mushy PCM exists in the central portion of the 
canister volume. Radial temperature gradients through this 
liquid zone are essentially zero giving rise to 99 percent 
side wall fractions through the 71 minute point in the 
cycle. Thereafter, the PCM freeze front advances radially 
outward beyond (r^ + r c )/2 establishing larger solid PCM 
temperature gradients in response to cooling fluid heat 
extraction at r^ This forces side wall fractions back into 
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the 30 to 40 percent range for the remaining portion of the 
cycle . 

5.3. 1.5 Limiting Effects of a Void 

The effect of a void on transient canister wall 
temperatures has not been quantified with the two- 
dimensional analyses described in this section. However, to 
derive a preliminary estimate of how a void could 
potentially increase predicted wall temperatures, the 
extreme case of a canister filled with PCM of zero thermal 
conductivity was analyzed. For this case, peak canister 
wall temperatures run between 20 K and 135 K higher than 
what is predicted for the canister with finite conductivity 
PCM. The likely peak wall temperature for a canister 
containing PCM with a void will be somewhere between the 
predictions of these two cases. Results from this analysis 
provide an upper limit of canister wall temperatures which 
are useful in developing the two-dimensional void heat 
transfer model. Canister thermal performance predictions 
with a void are discussed in the next section. 

5.3.2 Canister With Void Model 

5. 3. 2.1 Temperature and PCM Phase Distributions 


Figure 5.14 illustrates canister temperature contour and 
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(c) Time=69 . 80 min., MFL=0.5137, WF=0.1507, Nu=1.000 
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PCM phase maps at four times (24.28, 48.56, 69.80, and 81.94 
minutes) during a 91 minute orbital melt-freeze cycle. The 
PCM melting portion of the cycle occurs from time = 0 to "55 
minutes while PCM freezing occurs for the remaining portion 
of the cycle from time = “55 to "91 minutes. Initially, all 
PCM is solid at time = 0 minutes. During the PCM melting 
period, the large void thermal resistance forces a large 
percentage of canister heat transfer to occur via canister 
walls. This is illustrated in Figure 5.14a which shows high 
temperature gradients in the void region and isotherm 
normals (i.e., the direction of heat flow) generally aligned 
parallel to canister walls. Thus, energy absorbed in the 
canister outer wall diffuses around the void, down the side 
walls, and then into the PCM. As a consequence of this heat 
flow pattern, PCM melting occurs axially inward from both 
side walls. By time = 48.56 minutes when “90 percent of the 
PCM is liquid (see Figure 5.14b), heat transfer axially 
along the canister inner wall initiates PCM melting radially 
outward until all the PCM is liquified at time = “55 
minutes . 

PCM melting along the container walls is an extremely 
beneficial attribute of this TES canister design from a 
long-term structural integrity point of view. The two 
primary benefits include 1) structurally decoupling the 
solid PCM from the canister side walls and 2) providing a 
means through which expanding liquid PCM can flow into the 
void during the melting process. Decoupling the PCM from 
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the metal reduces canister thermal stresses created by 
differential thermal expansion between the PCM and metallic 
containment structure. PCM liquid flow paths to the void 
are highly desirable to preclude pressurizing entrapped 
liquid regions and the concomitant build-up of potentially 
large canister wall stresses. The importance of providing 
pressure relief flow paths is underscored by one 
experimental study involving a large volume of PCM (i.e./ a 
cube with 30 cm long sides), see Sparrow et al. (1978). For 
this experiment, the apparatus was constructed with a so- 
called vent heater installed in the PCM to guarantee a 
liquid flow path from the primary heated PCM region to the 
container void space. 

During PCM freezing, the heat of fusion energy liberated 
is transferred to the engine working fluid that cools the 
canister inner wall and to the canister outer wall where 
radiative heat loss to the receiver cavity occurs. Because 
the void acts as a thermal insulator, much of the heat loss 
from the liquid PCM to the canister outer wall occurs via 
conduction in canister side walls. As a consequence of this 
heat flow pattern, PCM freezing occurs along the canister 
inner wall and along the void surface. In addition, the 
maximum side wall temperature exists at about the radial 
midpoint (see Figure 5.14c). Near the end of the orbital 
cycle at time = 81.94 minutes, "90 percent of the PCM is 
frozen and the last remaining liquid PCM exists adjacent to 
void (see Figure 5.14d). Side wall radial temperature 
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profiles have been reestablished with temperature increasing 
in the positive radial direction. This is due to a small 
radiative heat input to the canister outer wall from the 
receiver cavity during the last "18 minutes of the orbital 
cycle (see Figure 3.3). 

5. 3. 2. 2 Void Heat Transfer 

Figure 5.15 illustrates void radial heat transfer, Q void , 
as a function of time during the orbital melt-freeze cycle. 
Q void is comprised of vapor conduction and surface-to-surface 
radiation components. By convention, these components are 
taken as positive if the resulting heat transfer is radially 
inward. During PCM melting, void heat transfer via 
radiation is about 3 times greater than that by vapor 
conduction and both components are positive and remain 
fairly constant. During PCM freezing, radiation is about 2 
times greater than vapor conduction and both components 
remain negative until "85 minutes into the cycle when all 
the PCM has frozen. The jump in the curves at "72 minutes 
is associated with the outer wall radiative heat flux 
boundary condition going from negative to positive. 

The fact that both components of void heat transfer 
remain negative from "72 to "85 minutes has interesting 
implications for the canister heat flow pattern. For this 
time period, relatively warm PCM transfers heat radially 
outward across the void to the canister outer wall where it 
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Figure 5*15* Void Radial Heat Transfer. 
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is then transferred back down the canister side wall and 
into the engine working fluid coolant (see Figure 5.14d). 

The conservative placement of a vapor-filled void volume 
adjacent to the canister outer wall generates the two 
negative effects of increased wall temperatures and 
increased wall temperature gradients (when compared to the 
case without void) during PCM melting. These effects 
increase canister thermal stresses and decrease canister 
design life predictions based on cumulative creep damage 
theory. Introduction of a void also reverses the sign of 
the side wall temperature gradient during the first half of 
PCM freezing period which does not occur in the case without 
a void. The resulting change in the canister side wall 
thermal stress distribution and the resulting impact (good, 
bad or indifferent) on canister life prediction is not 
easily determined without detailed structural analysis. 

A potentially beneficial effect of a void placed at the 
canister outer wall is a reduction in canister heat loss 
(more precisely, canister heat exchange with other canisters 
in the receiver cavity and heat loss out the cavity 
aperture) during PCM freezing due to the insulating quality 
of the void. The greatest canister heat loss occurs in the 
hottest canisters located near the aperture end of the 
receiver cavity (see Figure 2.2). These canisters are also 
located on the coolant tube near the inlet manifold and are 
thus cooled by relatively low temperature heat engine 
working fluid. Reduction of heat loss from these hottest 
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canisters permits greater heat transfer to the relatively 
cool working fluid near the inlet end of the coolant tube 
thereby decreasing the required working fluid heat transfer 
of canisters further down-stream. This, in turn, reduces 
the required temperature of down-stream canisters which must 
transfer heat to the highest temperature working fluid. 

5.3.3 Canister With Void and Free Convection Models 

5.3.3. 1 Temperature and PCM Phase Distributions 

Figure 5.16 illustrates the corresponding canister 
temperature contour and PCM phase maps for a 91 minute melt- 
freeze cycle which includes free convection in the liquid 
PCM. Figure 5.16a is identical to Figure 5.14a since the 
liquid PCM Rayleigh number (Ra) is below the critical Ra 
number and liquid PCM heat transfer is still controlled by 
conduction. Thus, no convective heat transfer enhancement 
takes place during the early part of the orbital cycle. At 
about 30 minutes into the cycle, the critical Ra number is 
exceeded and the Nu number begins steadily increasing from 
1.0 to a value of 4.5 at the end of the melting period ( 55 
minutes). Free convection in the liquid increases the rate 
of PCM melting and decreases canister temperature gradients 
as shown in Figure 5.16b where Nu = 3.047. During the PCM 
freezing portion of the cycle, liquid convective effects 
quickly die out and the Nu number falls back to 1.0 by 60 
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Figure 5.16. Canister Temperature Contour (Deg K) And 
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minutes into the cycle. Thus, for the majority of the PCM 
freezing period, liquid PCM heat transfer is again 
controlled by conduction. Hence, Figures 5.16c and 5.16d 
are nearly identical to Figures 5.14c and 5.14d with the 
exception that slightly more liquid PCM exists at the times 
of comparison for the case with free convection, i.e. the 
values of MFL (PCM mass fraction liquid) are slightly 
greater . 


5. 3. 3. 2 Effects of Free Convection 

Although relatively large Nu numbers exist during the 
PCM melting period, the length of time in which they occur 
is short. Thus over this short period of time, canister 
thermal performance in a 1-g environment, defined in terms 
of maximum wall temperature, canister temperature gradients, 
and PCM melting rate and/or PCM utilization, is not greatly 
different from that expected in a micro-g environment. For 
roughly 25 percent of canisters in the receiver cavity that 
contain high quality two-phase PCM (i.e., 0.5 < MFL < 1.0 ), 
a small reduction in maximum wall temperature and a small 
increase in PCM utilization would be expected during 1-g 
operation as a consequence of liquid PCM convection. For 
roughly 50 percent of the canisters that contain completely 
liquified PCM, approximately the same cyclic temperature 
range during 1-g and micro-g operation would be experienced. 
However, these canisters would experience a slightly lower 
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M time-at-maximum-temperature" history. For the remaining 25 
percent of canisters containing low quality two-phase PCM 
(i.e., 0.0 < MFL < 0.5 ), thermal performance would be 
essentially unaffected. Therefore, within the scope of this 
analysis, results indicate that canister thermal performance 
during ground tests should not be significantly different 
than that expected on-orbit. 

5.3.3. 3 Free Convection Model Assumptions 

These results must be interpreted in light of the 
assumptions of 1) axisymmetry (which requires alignment of 
the gravity vector and the canister axis of symmetry), 2) 
constant radiative flux input conditions at the canister 
outer wall in each case considered, 3) a prescribed void 
shape and location, and 4) axial-averaged liquid PCM 
characteristic length and radial temperature difference. 
Assumption 1 restricts the validity of analytical results to 
a ground test configuration with the canister axis vertical 
and with circumferentially uniform outer wall heating. 
Assumption 2 limits the available latitude for direct case- 
to-case comparison of results since large differences in 
canister temperature predictions invalidates the assumption 
of identical outer wall heat flux conditions in each case. 
This holds true since the canister outer wall boundary 
condition is a function of the radiation environment in the 
receiver cavity as well as canister outer wall temperature. 
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Assumption 3 removes the possibility of different void 
shapes and locations that are likely to occur between 1-g 
canister tests and micro-g canister operation. However, the 
canister performance differences associated with void shape 
and more importantly, void location, are bounded by results 
from two cases considered herein, i.e. a canister without a 
void and a canister with a void conservatively located 
adjacent to the outer wall. This assumption is further 
discussed in section 5.3.4. 

5. 3. 3. 4 Free Convection Model With Local Nu Numbers 

The validity of assumption 4) was checked by 
incorporating a local Nu number calculation into the two- 
dimensional canister computer program. For this 
calculation, the liquid PCM region characteristic length and 
radial temperature difference were determined for each PCM 
axial element group along the canister length from which a 
" z-dependent " Nusselt number, Nu(z), was calculated. Thus, 
the possibility for greater liquid PCM conductivity 
enhancement is permitted in locations adjacent to canister 
side walls where PCM liquifies first over centrally located 
PCM which liquifies last. 

Figure 5.17 shows Nu(z) as a function of time for the 
first 4 out of 8 PCM axial element groups. Nu(z) for the 
last 4 element groups are not shown since near longitudinal 
symmetry exists. The local Nu number first exceeds 1.0 at 
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2=0.2921 cm which is the coordinate of the PCM element group 
adjacent to the canister side wall at z=0. Further into the 
TES charging period, local Nu numbers exceed 1.0 at 
progressively larger z values until the last Nu(z) exceeds 
1.0 for the centrally located PCM element group at z=1.1303 
cm. The local Nu numbers individually increase to a plateau 
of about 5.5. This value corresponds to conditions of 
maximum characteristic length and moderated temperature 
difference due to the presence of two-phase PCM. At about 
52 minutes into the TES charging period, the 4 Nu(z) curves 
coalesce which corresponds to the time when complete PCM 
liquef ication has occurred. Nu(z) values then increase to a 
maximum value of about 7.5 at the end of the TES charging 
period in response to higher outer wall temperatures. 

Figure 5.18 shows a comparison of the previously 
calculated overall Nu number and the arithmetic mean of 
local Nu numbers calculated for one PCM charge-discharge 
cycle. The mean value of Nu(z), at any given instant, is 
about 60 to 70 percent higher than the overall Nu number 
calculated. This increase is attributed to the smaller 
effective convective cell aspect ratio, i.e. the PCM liquid 
height-to-width ratio. Yet the resulting impact on canister 
thermal performance is minimal. The predicted temperature 
and PCM phase distributions are nearly identical to those in 
Figure 5.16 in which the overall Nu number calculation was 
used. The predicted maximum canister wall temperature at 
any given time is not more than 5 K lower and the PCM mass 
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fraction liquid increases by only 1 to 2 percent. 

Therefore, the overall Nu number calculation method, which 
uses an averaged liquid PCM characteristic length and radial 
temperature difference, appears to be a satisfactory 
approach to further simplify modeling of convection effects 
in the liquid PCM region whose geometry is axially 
dependent . 

5.3.4 Performance Comparison 

In this section, results from three canister analytical 
cases are compared: 1) without a void, 2) with a void, and 

3) with a void and liquid PCM free convection (axial- 
averaged) . Comparison of results from cases 1 and 2 is 
intended to isolate the impacts of a void on canister heat 
transfer performance. Comparison of results from cases 2 
and 3 is intended to show likely canister performance 
differences during ground testing (in 1-g) and flight 
operation (in micro-g) . 

5. 3. 4.1 Maximum Wall Temperatures 

Figure 5.19 illustrates the maximum canister wall 
temperature, T(r,z) ma)( , throughout the 91 minute orbital 
cycle for cases 1, 2, and 3. T ( r / Z )*ax occurs at the axial 
midpoint of the outer wall and does not vary in position for 
the majority of the cycle. The introduction of a void 
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increases T(r,z) |WX about 20 K during the PCM melting period 
(0 to “55 minutes) and decreases T(r,z) (|a)( about 10 K during 
the PCM freezing period (“55 to “91 minutes). These changes 
are associated with the large void thermal resistance which 
increases outer wall axial temperature gradients by “100 
percent and increases side wall radial temperature gradients 
by *80 percent during the PCM melting period. 

The inclusion of liquid PCM free convection reduces 
T ( r , z ) Bax 0 to 10 K for the time between “30 and “55 minutes 
and has essentially no impact during other time periods (see 
Figure 5.19). During the same period, canister wall 
temperature gradients are “20 percent lower and the MFL is 
about 2 percent greater than for the case without free 
convection. Note that a further temperature reduction of up 
to 5 K is predicted when a local Nu is incorporated into the 
analysis. It is also interesting to note that although 
convective effects moderate the increase in T(r,z) MX during 
PCM melting, total PCM liquef ication occurs earlier than for 
the case without free convection. Thus, sensible heating of 
the liquid PCM occurs and quickly increases T(r,z) MX to 
about the same value that is predicted for the case without 
free convection. 

5. 3. 4. 2 Side Wall Heat Transfer Fractions 

Figure 5.20 illustrates the fraction of total canister 
radial heat transfer (comprised of void conduction, void 
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5.20. Radial Heat Transfer Side Wall Fraction. 
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radiation, and canister wall conduction) that occurs by 
conduction in the two canister side walls. This "side wall 
fraction" was evaluated at the radial location corresponding 
to the inner surface of the outer canister wall, R 0 '. The 
side wall fraction is a quantitative measure of how 
effectively canister walls redistribute energy absorbed in 
the outer wall and also provides a qualitative indication of 
wall temperature sensitivity to void location. 

The side wall fraction is relatively constant at *70 
percent for the case with a void (see Figure 5.20, "void" 
curve). Large perturbations in side wall fraction occur for 
brief periods of time at "55, *72, and "85 minutes. The 
three perturbations are associated with two step changes in 
the outer wall heat flux boundary condition and the point of 
complete PCM solidification, respectively. For the case 
without a void, the side wall fraction is considerably 
lower, i.e. 30 to 50 percent, due to the relatively high 
thermal conductance of the PCM when compared to the void. 

The same three perturbations exist as in the case with the 
void in addition to another perturbation at *62 minutes. 

This last perturbation is associated with a *7 minute period 
(from *55 to *62 minutes) during which a thin layer of PCM 
freezes along the outer wall. Once frozen, relatively large 
temperature gradients exist within the PCM layer which 
reduces the side wall fraction to only a couple percent at 
62 + minutes. 
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5. 3.4.3 Relationship Between Void Characteristics, Side 

Wall Fractions, and Wall Temperatures 

The prescribed void shape and location were chosen to 
permit a relatively straight-forward numerical analysis and 
to generate conservative predictions of canister wall 
temperatures. To confirm that void heat transfer and void 
placement are conservative, consider the relationship 
between maximum canister wall temperature and the radial 
heat transfer side wall fraction. As shown in Figures 5.19 
and 5.20, the canister wall radial heat transfer 
contribution can double (i.e., side wall fraction increases 
by "2 times) with only a 20 K increase in maximum wall 
temperature during PCM melting. Since the side wall 
fraction, by definition, is bounded by a 100 percent value, 
further increases in canister wall heat transfer are limited 
to “40 percent. Thus, regardless of the nature of void heat 
transfer, the increase in maximum canister wall temperature 
from introduction of a void is bounded by 28 K or a 
T ( r , 2 ) max of 1105 K for these case runs. This T(r,z) max value 
is just slightly higher than that predicted for the canister 
with void case (see Figure 5.19). Additionally, for any 
other viable void shape and location, PCM would be placed in 
contact with the outer wall thereby lowering the side wall 
fraction (see Figure 5.20, M no void" curve) and hence, 
lowering the maximum wall temperature. Therefore, 
differences in void shape and location between canister 1-g 
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ground tests and micro-g operation lead to relatively small 
changes in predicted canister thermal performance which are 
essentially bounded by the void cases considered herein. 

Results from two-dimensional canister analyses show that 
the largest changes in T(r,z) due to free convection and 
the introduction of a void are “-15 K and ~+20 K, 
respectively. These changes are an order-of -magnitude lower 
than the changes predicted from one-dimensional analyses in 
section 5.2.2 and the wall temperature increases from a void 
are 100 K lower than the preliminary estimates of section 
5.3.1. Therefore, this confirms the assertions from 
sections 5.2.2 and 5.3.1 that the effects of free convection 
and a void on canister thermal performance are much less 
pronounced in two-dimensional analyses than in a one- 
dimensional analyses. This result is attributed to the 
large heat transfer contribution of conduction within 


canister side walls. 



CHAPTER VI 


CONCLUSIONS 

6.1 One-Dimensional Analyses 

1. The present numerical solution methods accurately 
predict PCM canister thermal performance with and 
without the inclusion of a void (on the basis of 
comparisons with exact solutions). 

2. For a given thermal energy storage requirement, a 
constant charge/discharge boundary temperature results 
in lower boundary temperature variations from T m than 
with a constant flux boundary condition. 

3. Modes of void heat transfer must be selected to be 
consistent with the void size, void thermophysical 
properties, and the anticipated thermal environment. 
For the problems considered herein, void heat transfer 
from canister wall to PCM is best analyzed as a 
combined conduction-radiation process . 
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4. The presence of a one-dimensional void reduces material 
phase change rate under constant temperature boundary 
conditions by a factor of 4 to 5 and increases canister 
wall temperatures under constant flux input conditions 
by 50 to 200 K. 

5. The presence of liquid PCM free convection lowers 
canister wall temperatures by as much as 150 K and 
slightly enhances the PCM melting process during the 
heat input portion of the charge/discharge cycle. Free 
convection effects are essentially nonexistent during 
the TES discharge period. 

6. Ground-based receiver and canister performance in 1-g 
will be improved over on-orbit, micro-g performance. 
Therefore, ground test configurations which minimize 
free convective effects are suggested. The recommended 
test configuration places the canister and receiver 
axes of symmetry parallel to the gravity vector and the 
direction of heat input normal to the gravity vector. 
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6.2 Two-Dimensional A nalyses 


1. Based on the PCM freezing pattern from analyses without 
a void, a central core void location completely 
surrounded by solid PCM is anticipated at the end of 
the discharge period if PCM density change is taken 
into account. 


2. Canister walls effectively redistribute energy absorbed 
in the outer wall to the PCM and cooling fluid. 

Between 30 and 70 percent of the total canister radial 
heat transfer occurs by conduction within the canister 
side walls. 

3. Introduction of a void at the canister outer diameter, 
a) increases peak wall temperatures by 20 K, b) doubles 
canister wall temperature gradients during PCM melting, 
and c) transforms a predominately radial melting 
process into a predominately axial melting process . 


4. The void reduces radiative heat losses from the hottest 
canisters within the receiver cavity during the eclipse 
portion of the orbit. This beneficial effect improves 
the overall heat transfer to the engine working fluid 
and lowers canister temperatures during this period. 
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5. Void placement at the canister outer diameter is 
conservative and canister performance with any other 
viable void distribution is bounded by the performance 
predicted for the cases with and without a void 
considered herein. 

6. The presence of free convection in the liquid PCM: a) 

lowers peak wall temperatures by 0 to 15 K during the 
second half of the PCM melting period, b) lowers 
canister wall temperature gradients by 20 percent, c) 
increases the mass fraction of liquid PCM by ~2 
percent, and d) does not significantly alter the PCM 
melting/freezing process. 

7. An axial-averaged Nu number approach for modeling 
liquid PCM free convection heat transfer effects is an 
acceptable engineering approximation to a more 
detailed, local Nu number approach. 

8. One-dimensional analyses, that neglect canister side 
wall conduction, greatly over predict canister wall 
temperature changes associated with a void and free 
convection in the liquid PCM. 
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9. The differences in predicted canister performance 
between ground-based tests, in 1-g, and flight 
operation, in micro-g, are small on the basis of free 
convection heat transfer enhancement and void 
distribution effects. 

Future Work 


There are several options for extending the TES 
canister analytical work documented herein. Three 
particularly important options include: 1) further 

verification of analytical predictions, 2) extension of 
analytical methods to a three-dimensional canister geometry, 
and 3) parametric sensitivity and optimization studies. 

Analysis verification is a planned activity which will 
guantify the validity of the many engineering assumptions 
used to analyze canister heat transfer. This activity will 
utilize ground-based test data currently being generated at 
NASA Lewis Research Center, flight test data to be available 
in mid-1993 (see Namkoong (1989)), numerical predictions 
from the computer program developed by Wilson and Flanery 
(1988), and other numerical solutions available in the 
literature for idealized phase change problems. 

Extension to a three-dimensional geometry is necessary 
to predict non-axisymmetric canister thermal performance. 

The asymmetries are associated with: a) a non-uniform outer 

wall flux boundary condition (i.e., one half of the canister 
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circumference is exposed to direct concentrated solar energy 
and thermal radiation exchange with the receiver cavity 
while the other half of the circumference is exposed to 
reflected and reradiated energy from the receiver internal 
cavity wall) and b) ground testing in an orientation that 
does not maintain the canister axis of symmetry aligned 
parallel with the gravity vector. 

Parametric sensitivity studies are necessary to 
quantify the presumably small influence of uncertainties on 
analytical predictions. Uncertainties associated with the 
following items require consideration: material properties, 

boundary conditions, and the effective thermal conductance 
of the braze joint that attaches the canisters to cooling 
fluid tubes. Lastly, optimization studies could be 
performed to further refine the canister design. Variables 
on which to optimize include canister mass, PCM utilization, 
wall temperatures, and wall temperature gradients (i.e., 
thermal stresses). 
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The finite-difference energy balance equations for six 
canister model elements is given in Appendix Al. The six 
elements are located in the following positions: 1) in the 

central PCM ( l<i<ivl, l<j<jj ), 2) in the PCM at the PCM- 
void interface ( i=ivl, l<j<jj ), 3) in the canister outer 
wall ( i=ii, 1 < j < j j ), 4) in the canister side wall adjacent 
to the void ( ivl<i<ii, j-1 ), 5) in the canister side wall 
adjacent to PCM ( l<i<ivl, j=jj )/ and 6) in the canister 
inner wall ( i=l, l<j<jj ). The element numbering system is 
shown in Figure 4.1 and the symbols used in the following 
equations are defined in the nomenclature and/or in Appendix 
A2 , section A2 . 2 . 


Element 1 

E ( ,r =E M n 

Element 2 

+c3 ., 1 * k 5Pi,i"*t T «,i*' n - T *,i r ' 1 

-econ. n -qradj_ 1 n *At . 


PRECEDING PAGE BLANK NOT FILMED 


(A.l) 


(A. 2) 
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Element 3 


+c4 H,i* k 3'"ii,i"*[ T ii,j-i r '- T (,,i n ] 
+econ^ n +q*C10-qrad i+ . -_ 3 n *At . 

Element 4 

E i.r ,=E <.r +ci >,j* ki Pf, i "*[ T i.i, ) , '- T ,,j") 

+c 3 i,j* k jp>,r*t T <.i.i , ’- T t,i n ) 

-qrad^j j.j n *At . 


Element 5 


E l,)"*’ =E l.i" +C1 t,J* ki P.,i n *[T„ 1 , i n -T, j"] 

+C2 u* kim u n ’[ T i-i.) r '-' T .,,"I 

Element 6 

+ci i.)* ki P(,)"*t T »i.i , '- T -,j n i 

+C3 .,j* k iP.,)"*t T .,,*i n - T i,i"] 

+C6 *[T, "-T. "] . 

* J ■ t J 

In the above equations, the coefficient terms and the 
conductor terms are defined by: 

cl i, J " 2 * 7I * Az j* At / ln ( r i + i/ r i] / 


C2 iJ = 2 * rT * AZ j* At / ln t r i/ r i-1 1 i 


C3, j =4*n*r i *Ar i *At/[A2 jt1 +Az j ] , 


(A. 9) 
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C4, i =4*rr*r>Ar,*At/[Az 1 +Az. .] , 

’/J 11 J J 


ki P.,)" = W* k i,j"* I Ar «.1 +Ar i 1 /t k ) .1, ) "* to . +k 1,i‘* 4r «1 1 ’ 

kim l<t »-k l . 1 , j "*k t ,,**(Ar t . 1+ 4r 1 ] / [ k t . 1 , J "*4r,+k, <| «*Ar 1 . 1 ] , 


k iP.,j" =k 1 ,i*i n * k i,i"*[ Az i*i +Az i 1/[k ^i-i"* Az ) +k '-i"* Az i* 11 ' 


k 3”>.,i , ’* k .,i-<"* k <.i"* 1 Az i-i +Az i 1 -'l k i,i-i n * Az i +k <.i"* Az i-' 1 ' 


(A. 10) 
(A. 11) 
(A. 12 ) 
(A. 13) 
(A. 14) 


and where, 


C8=2*rc*ril*U*Az j *At 


(A. 15) 


C10 = 2*rc*ro2*Az ;| *At . 


(A. 16) 



Appendix A2 . FORTRAN Program Description and Listing 


FORTRAN 77 computer programs were written to analyze 
one- and two-dimensional containment canister heat transfer. 
Appendix A2 contains the listing and variable definitions 
for the two-dimensional computer program called NUCAM2DV 
which is an abbreviation for Numerical Canister Model: Two- 
dimensional with Void. A block diagram showing the main 
program and subroutines of NUCAM2DV is given in Figure A.l. 
The main program, MAIN . reads input data and executes a 
"time-marching" analysis via calls to various subroutines 
which are briefly described in the following paragraph. 

The subroutine INIT is called once at the beginning of 
program execution to read initial canister temperatures and 
initialize all program variables. Once in the time loop, 
calls are made at each time step to five subroutines. 

VOIDCON calculates steady state void vapor conduction heat 
transfer and temperature distributions. VOIDRAD calculates 
void surface element net radiation heat loss terms assuming 
that the void is a diffuse, gray, and opaque enclosure with 
a non-participating void vapor. SHELL performs an energy 
balance on the containment canister shell and determines 
canister wall temperature distributions. FLUID calculates 
he/xe cooling fluid pseudo-steady state temperature 
distributions and total heat transfer to the cooling fluid. 
SALT performs an energy balance on the LiF-CaF 2 salt PCM and 
determines PCM temperatures, phase distributions, and up- 
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(9 FUNCTION 
SUBROUTINES) 



Figure A.l. NUCAM2DV Computer Program Block Diagram 
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dated material properties. The last subroutine, UPDATE , is 
called by MAIN at two different time increments during 
program execution (typically, 1.0 min. and 3.035 min.). 
UPDATE , in turn, makes calls to five subroutines which are 
briefly described below. 

VOIDFF calculates void surface element view factors 
(using closed-form solutions and view factor algebra) and 
inverts the "view factor - emittance" matrix. From VOIDFF . 
"calls” to nine different function subroutines are made. 
Function subroutines are needed to evaluate view factor 
algebra equations which employ multiple calls from the same 
line of code. CONV calculates liquid PCM Rayleigh numbers, 
Nusselt numbers, and enhanced conductivity values. ENERGY 
evaluates the global canister energy balance. WALLFRAC 
calculates side wall radial heat transfer fractions near the 
inner, mean, and outer canister radial locations. The last 
subroutine, OUTPUT , creates files for numerical data output. 



on on 
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A2.1 Two-Dimensional Analysis Program Listing 
C 

C 2-D (R,Z), TRANSIENT TES CANISTER MODEL: NUCAM2DV 
C 

C INCLUDES VOID AND FREE CONVECTION AS FUNCTION OF Z 
C 

C THOMAS W. KERSLAKE NASA-LEWIS RESEARCH CENTER CLEVELAND, OHIO 
C SOLAR DYNAMIC POWER SYSTEM BRANCH 

C 

C ADVISER DR MOUNIR B IBRAHIM CSU DEPT. OF MECH. ENG 
C 

C START DATE: 20 NOV 89 REV DATE23 JULY 90 
C 

C SIMPLE EXPLICIT NUMERICAL FORMULATION / ENTHALPY METHOD 
C 

C UNITS : LENGTH CENTIMETER 

C MASS GRAM 

C ENERGY JOULE 

C TIME SECOND 

C TEMPERATURE DEG K 

C 

CHARACTERS PHASE(40,20),NC 

REAL MDOT 1 KW,KS 1 KL I KLE,MT,MS 1 ML,K(40 1 20),KIP(40 I 20) 1 
&K1M(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU I NUL,KF 

COMMON DT.DRR.DZZjDELWO.DELWI.DELWS.Rll.WP.BBINVflS.lBJ.QRADllS), 
&RI2,R01,R02,Q,TFI,MD0T,CPF,KW,CPW,RH0W,KS,CPS,RH0S,FF(18,18), 
8iKL,CPL,RHOL,TM,HM,KLE,Z2,Z3,Z4,PI,VOLT,Ql,Q2,VOL(40,20),RCML(20), 
&MT,MS,VOLS,ML,VOLL,ESHSUM,EFSUM,RA,NU,CHARL,C7,C12,ECON(20), 
&TV,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20),AREAOLD(3), 

AYF (40,20) , DR(40) ,DZ(20) ,Z(20) , R(40) ,EO(40,20) , E(40,20),AREA(30), 

4£TFO(20),TF(20),C1(40,20),C2(40,20),C3(40,20),C4(40,20),RCMLO(20), 

&TFOUT,EE(40,20) 1 TIMEPR1,TIMEPR2,TIMEP,TIMEF,TFP,TFF,RV,RM, 

&TIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL,EEL(40,20), 

fcEELO(40,20),EEL02(40,20),PR,EML(40,20), 

&MFL,QF, II, IV1,IV2,JJ,N, PHASE, NC, 

&ZKLE (1 0) ,ZNU (1 0) ,ZRA(1 0) ,ZCHARL(1 0) ,ZNUAVG 
C 

INPUT DATA 

READ (1,1000) DT,NN,DRR, II, DZZ.JJ.DELWO, DELWI.DELWS.RIl, 

St NC,MT,MDOT,PRF,MUF,CPF,KF,RHOF,KW,CPW,RHOW,EW, 

St KS, CPS, RHOS,ES,KL,CPL,RHOL,TM,HM,NUL, BETA, RAC, 

St G.TIMEPR1 .TIMEPR2, 

St TIME1,TIME2,TIME3,TIME4,TF1,TF2,TF3,TF4, 

A Q1.Q2.Q3 

CALCULATE PROGRAM CONSTANTS 
C 
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AS=KS/(RHOS*CPS) 

AL=KL/(RHOL*CPL) 

PR=NUL/AL 

GAMMAO = 1 . DOO-RHOL/RHOS 
KLE=KL 

RI2=RI1+DELWI 
ROl = RI2+ DRR*(ll-2) 

R02=R01+DELW0 

Z2=DELWS 

Z3=DELWS+DZZ*(JJ-2) 

Z4=Z3+DELWS 

JJ2=JJ/2 

TIMEP=TIME1 

TIMEF=TIME2 

TFP=TF1 

TFF=TF2 

PI =2 ODOO‘DARSIN(1.0DOO) 

V=MDOT/(RHOF*PI*Rll**2) 

RE=2 ODOO*MDOT/(PI*MUF*RI1) 

C 

C CONSTANT PROPERTY, FULLY DEVELOPED FLOW CORRELATION (KAYS 1966) 
C 

H=KF*0 022D00*RE**0 8D00*PRF**0 6D00/(2.ODOO*Rll) 

U = 1 D00/(1 DOO/H + DLOG((Rll + DELWI/2 D00)/RI1)*RI1/KW) 

DUM=PI 

TFI=TF1 

Q=Q1 

C 

C INITIALIZE ARRAYS AND VARIABLES TO AN "ALL SOLID SALT' STATE 
C FOR WHICH ALL CANISTER TEMPERATURES ARE < TM AND DEFINE 
C GEOMETRY ARRAYS 
C 

CALL INIT(WF,VOLV) 

IF (RV.LT.ROl) CALL VOIDFF(ES,EW) 

C7=G*BETA/(AL*NUL) 

C8=2 0D00*PI*RI1»U*DT*DZZ 

C9=2 0D00*PI*RI1*U*DT*DELWS 

C10=2 0D00*PI*RO2*DT*DZZ 

Cll=2.0D00*PI*RO2*DT*DELWS 

C12=2 0D00*PI*RO2*Z4 

C13=U*PI*RI1*DELWS 

C14=MDOT*CPF 

C15=C13+C14 

C16=U*PI*RI1*DZZ 

C17=C16+C14 

MS=MT 

VOLS = MS/RHOS 
QSHELL=Q*C12 
TTUBE = 0 0D00 
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DO 11 J = 1,JJ 
TTUBE=TTUBE+TO(l,J) 

11 CONTINUE 

TTUBE=TTUBE/JJ 

QF=U*2 D00*PI*RI1*Z4*(TTUBE-TFI) 

QSALT=QSHELL-QF 

TV=TO(IVl,JJ2) 

NU = 1 ODOO 
ZNUAVG = 1 ODOO 
QVOIDCON = 0 ODOO 
QVOIDRAD=0 ODOO 
QVOIDTOT=O.ODOO 
C 

C ECHO INPUT AND INITIALIZED VARIABLES 
C 

WRITE(7,1000) DT,NN,DRR,II,DZZ,JJ,DELW0,DELWI,DELWS,RI1, 

St NC,MT,MDOT,PRF,MUF,CPF,KF,RHOF,KW,CPW,RHOW,EW, 

Si KS, CPS, RHOS,ES,KL,CPL,RHOL,TM,HM,NUL, BETA, RAC, G.TIMEPRl, 

Si TIMEPR2,TIME1,TIME2,TIME3,TIME4,TF1,TF2,TF3,TF4, 

& Q1,Q2,Q3 

WRITE (7,2000) AS,AL,PI,GAMMA0,PR,RI2,RV,R01,R02,Z2,Z3,Z4, 

St VOLT,VOLS,VOLL,VOLV,MT,MS,ML,V,RE,H,U 

WRITE (8,2100) DT,NN,DRR,II,DZZ,JJ,DELW0,DELWI,DELWS,RI1,R02,Z4, 

Si MT,MD0T,H,NC,Q1,Q2,Q3 
WRITE (8,2200) TIMEM,MFL,T(II,JJ2),WF, 

4£TFI,TFOUT,QF,EBAL,EBALP,QVOIDCON,QVOIDRAD,QVOIDTOT 

WRITE (11,2410) 

WRITE (11,2420) TIMEM,WALLFR1,WALLFR2,WALLFR3 
IF(NC EQ ’ON ’) WRITE(9,2300) DT.NN.DRR.II.DZZ.JJ.DELWO, 

Si DELWI,DELWS,MT,TFI,MD0T 1 H,NC,RAC,Q1,Q2,Q3 
IF(NC EQ 'ON ’) WRITE (9,2400) TIMEM,CHARL,TSHELL,RA,NU 
IF(NC EQ 'ON ’) WRITE (17,3600) 

IF(NC.EQ ’ON ’) WRITE (18,3605) 

IF(NC EQ ’ON ') WRITE (17,3650) TIMEM,ZNUAVG,(ZNU(J),J=2,JJ-1), 

Si (ZCHARL(J),J=2,JJ— 1) 

IF(NC EQ ’ON ’) WRITE (18,3655) TIMEM,ZRAAVG,(ZRA(J),J=2,JJ-1) 

N=2 

CALL 0UTPUT(TIMEM,N1,N2,WALLFR1,WALLFR2,WALLFR3,EBALP,WF) 

START SIMULATION 

N1 = (TIMEPR1/DT) + 1 0 
N2 = (TIMEPR2/DT) + 1.0 
N4=2 

DO 20 N-2,NN 

C STEADY STATE VOID CONDUCTION HEAT TRANSFER IN RADIAL DIRECTION 


CALL VOIDCON(DUM) 



150 


C 

C DIFFUSE, OPAQUE, AND GRAY VOID RADIATION HEAT TRANSFER 
C 

CALL VOIDRAD(DUM) 

C 

C ENERGY BALANCE ON SHELL (CONTAINMENT CANISTER) ELEMENTS 
C 

CALL SHELL(C8,C9,C10,C11) 

C 

C ENERGY BALANCE ON HE/XE GASEOUS COOLING FLUID ELEMENTS 
C 

CALL FLUID(C13,C14,C15,C16,C17) 

C 

C ENERGY BALANCE ON EUTECTIC LIF-CAF2 SALT PCM 
C 

CALL SALT(DUM) 

C 

C "ADVANCE" TIME-DEPENDENT VARIABLES TO THE "N+l" TIME STEP 
C 

DO 25 J = 1,JJ 
TFO(J)=TF(J) 

DO 30 1 = 1,11 
TO(l,J)=T(l,J) 

EO(l,J) = E(l,J) 

30 CONTINUE 
25 CONTINUE 
C 

C UP-DATE VOID EFFECTS, FREE CONVECTION EFFECTS, ENERGY BALANCE 
C TALLY, WALL FRACTION CALCULATION, AND OUTPUT RESULTS 
C 

IF(N EQ N1 OR N EQ N2) CALL UPDATE(N1,N2,WF) 

C 

C OUTPUT TEMPERATURE DATA FOR VIDEO ANIMATION 
C 

IF (N NE N4) GO TO 20 
TIME=(N-1)*DT 

WRITE (16,2450) TIME,MFL,WF,ZNUAVG,(Z(J),J=1,JJ) 

WRITE (16,2050) 

DO 855 1=1,11 
111 = 11 + 1-1 

WRITE(16,2500) R(II1),(T(II1,J),J=1,JJ) 

855 CONTINUE 

N4=N4+6 0D00/DT 
20 CONTINUE 
C 

1000 FORMAT(F9 7,1X,I6,2(1X,F7.5,1X,I3),4(1X,F7 5)/ 

& A4.1X.F6 3,1X,F6 3,1X,F7.5/ 

* 4(D9 3,1X)/4(F7 5,1X)/4(F7 5,1X)/3(F7 5,1X),F61,1X, 

L F5 1,1X,F6 4,1X,D8.2,1X,F7.1,1X,F5.1/F6.1,1X,F6.1/ 
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Sc 4(F6 2,1X)/4(F6 1,1X)/3(F6 3, IX)) 

2000 FORMAT(6(5(D12 5,1X)/)) 

2100 FORMAT(’DT=’,F9.7,’ NN = ’,I6,’ DR=\F7.5,' 11=’, 13,’ DZ=',F7 5, 

Sc ’ JJ = ’ 1 I3/ , DELW0=’ I F7.5,’ DELWI = ’,F7.5,’ DELWS= , ,F7.5 1 ’ Rl=’, 

&F7.5,* RO=’,F7 5/’L= ',F7 5/ MT=’, 

*F6.3,’ MDOT=’,F6 3,’ H=’,F75,’ NC=’,A4/ 

&*Q1 = ’,F6 3,’ Q2=’ 1 F6.3 1 * Q3=’,F6.3/ 

Sc /’ TIME'.TlO.’MFL'/rie.’TMAX’/ra,’ WF ’,T33, 

&’TFr,T40,’TFO’,T49,’QF’ 1 T56,’EBAL’,T63,'%EBAL’ 1 

4£2X j 'QVOIDCON’,2X,'QVOIDRAD’ j 2X,’QVOIDTOT7/) 

2200 FORMA^Fe^^.FB.S.lX.Fe i^.Fe^^ZX^e.lj.lX.FZZ, 

Sc 2X,F7.0,2X 1 F5.2 J 2X 1 3(F8.3,2X)) 

2300 FORMAT( , DT=’ l F9.7 l ’ NN = ’,I6,’ DR=’ I F7.5,’ 11=’, 13/ DZ=’,F7.5, 

Sc ' JJ = ’,I3/ 1 DELW0=’,F7.5 1 ’ DELWI = ’,F7 5,’ DELWS=’ 1 F7.5 1 ’ MT =’ 

4£,F6.3/’ TFI(0 MIN) = , 1 F6.1,’ MDOT=*,F6 3 1 ’ H=’,F7.5 1 ’ NC=’,A4, 

A1X,’ RAC=',F71/ , Q1 = , ,F6.3 > ’ Q2= , ,F63,' Q3=',F63// 

Sc ' time’,t8,’charl’,ti4,’ tv , i t26,’Ra , ,T36 i *nuv) 

2400 F0RMAT(F6.2,1X,F5 3,2X,F6 1,2X,F9.0,2X,F6.3) 

2410 FORMAT(’ WALL HEAT TRANSFER FRACTIONS, %'//’ TIME ’,5X, 
&’WALLFR1’,3X,’WALLFR2’ I 3X,’WALLFR3’/) 

2420 FORMAT(F7 3,3(3X,F7 3)) 

2050 FORMAT(’R, CM’) 

2450 FORMAT(/’ TIME = ’,F6 1,’ MFL=\F6.4,’ WF = , J F6.4 I ’ NU=’,F63 
Sc //’ Z, CM = ’,20(2X,F6 4)) 

2500 FORMAT(F6.4,6X,20(F6 1,2X)) 

3600 FORMAT(’ TIME’.TS.’ZNUAVG’.TSO/NUlZJ’.TSO.’CHARLp)’) 

3605 FORMAT(’ TIME’.Tll/ZRAAVG’.TSO.’RAfZ)’) 

3650 FORMAT^e^.TX.FSJ.TX.SlFS S.lXj.lX.SlFS ^lX)) 

3655 FORMAT(F6 2,2X,F9 0,2X,8(F9.0,1X)) 

STOP 

END 

C 

SUBROUTINE INIT(WF,VOLV) 

CHARACTERS PHASE(40,20),NC 

REAL MDOT I KW I KS,KL,KLE 1 MT,MS,ML,K(40 1 20),K1P(40,20) ) 
&KIM(40,20) 1 KJP(40,20),KJM(40,20),MFL,MUF 1 NU,NUL J KF 

COMMON DT.DRR.DZZ.DELWO.DELWI.DELWS.RII.KIP.BBINVIIS.ISJ.QRADIIS), 
<iRI2,R01,R02,Q,TFI,MD0T,CPF,KW 1 CPW 1 RH0W,KS I CPS,RH0S,FF(18,18) l 
4 £ KL,CPL i RHOL i TM 1 HM,KLE 1 Z2 1 Z3 1 Z4,PI,VOLT,Q1 > Q2 1 VOL(40 ( 20) ( RCML(20) j 
4£MT ) MS,VOLS,ML,VOLL 1 ESHSUM,EFSUM 1 RA,NU,CHARL j C7 j C 12,ECON(20) 1 
1TV,TO(40,20),T(40,20),K 1 RHO(40,20) I XF(40 ) 20),TFIN(20) I AREAOLD(3), 
8 £ YF(40,20),DR(40),DZ(20) 1 Z(20),R(40),EO(40,20),E(40 1 20) 1 AREA(30), 
tTFO(20),TF(20),Cl(40,20),C2(40 J 20) 1 C3(40 l 20),C4(40 J 20) 1 RCMLO(20), 
1TFOUT,EE(40 1 20),TIMEPR1,TIMEPR2,TIMEP 1 TIMEF I TFP,TFF ) RV J RM, 
&TIME2,TIME3,TIME4,TF2,TF3,TF4 1 EINIT,ECAN J EBAL,EEL(40,20) J 
4'EELO(40,20),EEL02(40 1 20) J PR,EML(40 I 20) 1 
&MFL,QF,II,IV1 ) IV2,JJ 1 N 1 PHASE,NC 1 
42KLE(10),ZNU(10),ZRA(10),ZCHARL(10),ZNUAVG 
C 



n n 


152 


C READ IN INITIAL ELEMENT TEMPERATURES AND DETERMINE INITIAL ELEMENT 
DIMENSIONS, MATERIAL PROPERTIES, AND GEOMETRY COEFFICIENT MATRICES 

DO 41 1=1,11 
111 = 11 + 1-1 

READ (2,1100) (TO(lll,J),J=l,JJ) 

41 CONTINUE 
DO 25 1=2,11-1 
DR(I)=DRR 

R(l) = RI2+ DRR/2. D00+ (l-2)*DRR 
DO 30 J=2, JJ— 1 
ZKLE(J)=KL 
ZNU(J)=1 0D00 
K(I,J)=KS 
RHO(l,J)=RHOS 
PHASE(l,J) = ’SOL’ 

XF(l,J)=0 0D00 
YF(l,J)=0 0D00 
DZ(J)=DZZ 

VOL(I,J)=2,ODOO*PI*R(I)*DR(I)*DZ(J) 

EO(l,J) = (TO(l,J)-TM)*CPS*RHO(l,J)*VOL(l,J) 

30 CONTINUE 
25 CONTINUE 

R(l) =RI1 +DELWI/2.D00 
R(ll) = ROl + DELWO/2. D00 
DR(1)=DELWI 
DR(ll)=DELWO 
DO 35 J =2,JJ— 1 

Z(J) = DELWS+DZZ/2 D00+(J-2)*DZZ 

K(1,J)=KW 

K(II,J) = KW 

TFO(J)=TFI 

TFIN(J)=TFI 

VOL(l,J)=2 0D00*PI*R(1)*DR(1)*DZ(J) 

VOL(ll,J)=2 0D00*PI*R(ll)*DR(ll)*DZ(J) 

EO(l,J) = (TO(l,J)-TM)*CPW*RHOW*VOL(l,J) 

EO(ll,J) = (TO(ll,J)-TM)*CPW*RHOW*VOL(ll,J) 

35 CONTINUE 
DZ(1)=DELWS 
DZ(JJ)=DELWS 
Z(l)=DELWS/2 DOO 
Z(JJ) = DELWS/2 DOO + Z3 
DO 40 1=2,11-1 
K(I,1)=KW 
K(I,JJ) = KW 

VOL(l,l)=2 0D00*PI*R(I)*DR(I)*DZ(1) 

VOL(l,JJ)=2 0D00*PI*R(l)*DR(l)*DZ(JJ) 

EO(l,l) = (TO(l,l)-TM)*CPW*RHOW*VOL(l,l) 

EO(l,JJ) = (TO(l,JJ)-TM)*CPW*RHOW*VOL(l,JJ) 
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40 CONTINUE 
TFO(l)=TFI 
TFO(JJ)=TFI 
TFOUT=TFI 
TFIN(JJ)=TFI 

VOL(l,l)=2 ODOO*PI*R(1)*DR(1)*DZ(1) 

VOL(l,JJ)=2 0D00*PI*R(l)*DR(l)*DZ(JJ) 

VOL(ll,l)=2 0D00*PI*R(II)*DR(II)*DZ(1) 

VOL(ll,JJ)=2 0D00*PI*R(ll)*DR(ll)*DZ(JJ) 

EO(l,l) = (TO(l J l)-TM)*CPW*RHOW*VOL(l,l) 
EO(l,JJ) = (TO(l 1 JJ)-TM)*CPW*RHOW*VOL(l,JJ) 
EO(II 1 1) = (TO(II,1)-TM)*CPW*RHOW*VOL(II,1) 

EO(ll, JJ) = (TO(ll,JJ)-TM)*CPW*RHOW*VOL(ll 1 JJ) 
K(1,1) = KW 
K(1,JJ) = KW 
K(II,1) = KW 
K(II,JJ)=KW 
C 

C RECALCULATE VOID AFFECTED VARIABLES 
C 

IV1 = II-1 
IV2=IV1 

RV=DSQRT(1 ODOO*RI2*RI2+MT/(RHOS*PI*(Z3-Z2))) 
IF (RV EQ ROl) GO TO 44 
11 = 11-1 

IV1 = (RV-RI2) /DRR+1 0D00 
IV2=IV1 + 1 

DR(IVl) = RV-(R(IV1 -1) + DRR/2. D00) 

DR(IV2) = (R(IV2+2)-DRR/2.D00)-RV 
DR(ll) = DELWO 
R(IVl) = RV-DR(IVl)/2 0D00 
R(IV2) = RV+ DR(IV2)/2 0D00 
DO 37 J=2,JJ-1 
K(II,J) = KW 
DO 36 I = IV2.II — 1 

K(I,J)=4 6D-04*DSQRT(TO(I,J)/1000 0D00) 
PHASE(l,J)=’VOID’ 

EO(I,J)=0 0D00 
E(l,J)=0 0D00 

36 CONTINUE 

VOL(IV1,J)=2 0DOO*PI*R(IV1)*DR(IV1)»DZ(J) 

E0(IV1, J) = (T0(IV1 , J)-TM)*CPS*RHOS*VOL(IVl ,J) 

37 CONTINUE 
DO 38 I = IV1,II 

R(l) = R(l-l) + (DR(I) + DR(l-l))/2 ODOO 
VOL(l,l) =2 0D00*PI*R(I)*DR(I)*DZ(1) 

EO(l,l) = (TO(l,l)-TM)*CPW*RHOW*VOL(l,l) 
VOL(l,JJ)=2 0D00»PI*R(l)*DR(l)*DZ(JJ) 

EO(l,JJ) = (TO(l,JJ)-TM)*CPW*RHOW*VOL(l,JJ) 
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38 CONTINUE 

DO 39 J =2,JJ — 1 

VOL(ll,J)=2.0D00*PI*R(ll)*DR(ll)*DZ(J) 

EO(II.J) = (TO(ll, J)-TM)*CPW*RHOW*VOL(ll, J) 

39 CONTINUE 
R01=R(ll)-DR(ll)/2 0D00 
RO2=R(ll) + DR(ll)/2 0D00 

44 00 45 1=2,11-1 
DO 50 J=2,JJ-1 

C1(I,J)=2 DOO*PI*DZ(J)*DT/DLOG(1.0DOO*R(I+1)/R(I)) 
C2(I,J)=2 DOO*PI*DZ(J)*DT/DLOG(1,ODOO*R(I)/R(I— 1)) 
C3(I,J)=4D00*PI*R(I)*DR(I)*DT/(DZ(J + 1)+DZ(J)) 
C4(I,J)=4D00*PI*R(I)*DR(I)*DT/(DZ(J-1)+DZ(J)) 

50 CONTINUE 

45 CONTINUE 
DO 55 1=2,11-1 

Cl (1, 1) =2.DOO*PI*DZ(1)*DT/DLOG (1 ODOO*R(l + 1)/R(I)) 
C2(l,l)=2 DOO*PI*DZ(1)*DT/DLOG(1.0DOO*R(I)/R(I— 1)) 
C3(l,l) = 4 D00*PI*R(l)*DR(l)*DT/(DZ(2) + DZ(l)j 
C1(I,JJ)=2 DOO*PI*DZ(JJ)*DT/DLOG(1,ODOO*R(I+1)/R(I)) 
C2(I,JJ)=2 DOO*PI*DZ(JJ)*DT/DLOG(1 ODOO*R(l)/R(l-l)j 
C4(I,JJ)=4D00*PI*R(I)*DR(I)*DT/(DZ(JJ-1)+DZ(JJ)) 

55 CONTINUE 
DO 60 J =2, JJ— 1 

Cl (1 , J) =2 D00*PI*DZ(J)*DT /DLOG(l 0D00*R(2)/R(1)) 
C3(1,J)=4D00WR(1)*DR(1)*DT/(DZ(J+1)+DZ(J)) 
C4(1,J)=4D00*PI*R(1)*DR(1)*DT/(DZ(J-1)+DZ(J)) 
C2(II,J)=2 DOO*PI*DZ(J)*DT/DLOG(1 0D00*R(II)/R(II-1)) 
C3(II,J)=4D00*PI»R(II)*DR(II)*DT/(DZ(J + 1)+DZ(J)) 
C4(II,J)=4D00*PI*R(II)*DR(II)*DT/(DZ(J-1)+DZ(J)) 

60 CONTINUE 

C1(1,1)=2DOO*PI*DZ(1)»DT/DLOG(1,ODOO*R(2)/R(1)) 
C3(l,l) = 4 D00*PI*R(1)*DR(1)*DT/(DZ(2)+DZ(1)) 

C 1 (1 , J J) = 2 D00* PI*DZ(J J) * DT/DLOG (1 0D00*R(2)/R(1)) 
C4(l,JJ)=4D00*PI*R(l)*DR(l)*DT/(DZ(JJ-l)+DZ(JJ)j 
C2(ll,l)=2 DOO*PI*DZ(1)*DT/DLOG(1 0D00*R(II)/R(II-1)) 
C3(II,1)=4D00*PI*R(II)»DR(II)»DT/(DZ(2) + DZ(1)) 
C2(II,JJ)=2 D00*PI*DZ(JJ)*DT/DLOG(1.0D00*R(ll)/R(ll— 1)) 
C4(II,JJ)=4D00*PI*R(II)*DR(II)*DT/(DZ(JJ-1)+DZ(JJ)) 
V0LV=PI*(R01*R01-RV*RV)*(Z3-Z2) 
V0LT=PI*(R01**2-RI2**2)*(Z3-Z2) 

WF = VOLV/VOLT 
EINIT=0 0D00 
DO 65 1=1,11 
DO 70 J=1,JJ 
EINIT=EINIT +EO(l,J) 

70 CONTINUE 
65 CONTINUE 
ECAN = EINIT 
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1100 FORMAT(12X,20(F6 1 ,2X)) 

RETURN 

END 

SUBROUTINE VOIDCON(DUM) 

CHARACTERS PHASE(40,20),NC 

REAL MDOT,KW,KS,KL,KLE,MT,MS,ML,K(40,20),KIP(40,20), 
4dKIM(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU,NUL,KF 

COMMON DT.DRR^ZZ.DELWO.DELWI.DELWS.Rll.WP.BBINVtlS.lSj.QRAD^S), 
&RI2,R01,R02,Q,TFI,MD0T,CPF,KW,CPW J RH0W,KS,CPS,RH0S,FF(18,18), 
8iKL,CPL,RHOL 1 TM 1 HM,KLE J Z2 ) Z3,Z4,PI J VOLT,Ql,Q2,VOL(40 l 20) ) RCML(20), 
&MT,MS,VOLS,ML,VOLL,ESHSUM,EFSUM,RA,NU,CHARL 1 C7,C12 I ECON(20) J 
<n"V,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20),AREAOLD(3) J 

4tYF (40 , 20) , DR (40) , DZ(20) , Z(20) , R (40) , EO (40,20) , E (40,20) , AREA(30) , 

&TFO(20),TF(20),C1(40,20),C2(40,20),C3(40,20),C4(40 1 20) 1 RCMLO(20), 

&TFOUT,EE(40,20),TIMEPR1,TIMEPR2,TIMEP,TIMEF,TFP I TFF ) RV,RM, 

4tTIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL 1 EEL(40 I 20), 

fcEELO(40,20),EEL02(40,20),PR,EML(40,20), 

4tMFL,QF, II, IV1,IV2,JJ,N, PHASE, NC, 
&ZKLE(10),ZNU(10),ZRA(10),ZCHARL(10),ZNUAVG 
DO 10 J =2, JJ — 1 

C=(T0(II,J)-T0(IV1,J))/DL0G(R01/RV) 

D=T0(IV1,J)-DL0G(RV)*(T0(II,J)-T0(IV1,J))/DL0G(R01/RV) 

DO 20 l = IV2,ll-l 
T(I,J) =C*DLOG(R(l)) + D 
TO(l,J)=T(l,J) 

K(I,J)=4 6D -04*DSQRT(T(I,J)/1000 0D00) 

20 CONTINUE 
10 CONTINUE 
RETURN 
END 
C 

SUBROUTINE VOIDRAD(DUM) 

CHARACTERS PHASE(40,20),NC 

REAL MDOT,KW,KS,KL,KLE,MT,MS,ML,K(40,20),KIP(40 > 20) > TTO(18), 
AKIM(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU,NUL 1 KF,DD(18) 

COMMON DT.DRR.DZZ.DELWO.DELWI.DELWS.Rll.KIP.BBINVjlS.lSj.QRADflS), 
<tRI2,R01,R02,Q,TFI,MDOT,CPF,KW,CPW,RHOW,KS,CPS J RHOS ) FF(18,18), 
&KL,CPL,RHOL,TM,HM,KLE,Z2,Z3,Z4,PI,VOLT,Q1,Q2,VOL(40,20),RCML(20) 1 
&MT,MS,VOLS,ML,VOLL,ESHSUM,EFSUM,RA,NU,CHARL,C7,C12 1 ECON(20), 
4tTV,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20) ) AREAOLD(3) J 
4 £ YF(40,20),DR(40),DZ(20),Z(20),R(40),EO(40,20),E(40 ) 20) 1 AREA(30) J 
tTFO(20),TF(20),Cl(40,20),C2(40,20),C3(40,20),C4(40 1 20),RCMLO(20) l 
1TFOUT,EE(40,20),TIMEPR1,TIMEPR2,TIMEP,TIMEF 1 TFP J TFF 1 RV I RM, 
iTIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL,EEL(40 1 20) 1 
1EELO(40,20),EEL02(40,20),PR,EML(40,20), 

&MFL.QF, II, IV1,IV2,JJ,N, PHASE, NC, 

4tZKLE (1 0) ,ZNU (1 0) ,ZRA(1 0) ,ZCH ARL(IO) ,ZNU AVG 
SIGMA=5.67D-12 
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NRS=2*(JJ-2)+2 
DO 5 KK=l,JJ-2 
TTO(KK)=TO(IVl,KK+l) 

TTO(KK+ JJ-2) =TO(ll,KK+ 1) 

5 CONTINUE 
l=IV2 

JRANGE=2*(JJ-2) + 1 
TTO(JRANGE) =0 0D00 
TTO(JRANGE+ 1) =0 0D00 
DO 6 KK=JRANGE,NRS-1,2 
DO 7 I = IV2.II— 1 
TTO(KK) =TTO(KK) +TO(l,l) 

TTO(KK+ 1) =TTO(KK+ 1) +TO(l, JJ) 

7 CONTINUE 

TTO(KK) =TTO(KK)/(ll-IV2) 

TTO(KK+ 1) =TTO(KK+ 1)/(II-IV2) 

6 CONTINUE 
C 

C DEFINE "DD" MATRIX SUCH THAT {BB}{QRAD} = {DD} WHERE {BB> IS 
C THE 'VIEW FACTOR-EMITTANCE" MATRIX AND {DD} IS THE 
C "TEMPERATURE-TO-THE-FOURTH-POWER DIFFERENCE" MATRIX 
C 

DO 10 KK=1,NRS 
SUM=0 0D00 
DO 20 J = 1,NRS 

SUM=SUM+FF(KK,J)*(TTO(KK)*M-TTO(J)**4) 

20 CONTINUE 

DD(KK)=SIGMA*SUM 
10 CONTINUE 
C 

C DETERMINE ELEMENT NET HEAT LOSS TERMS 
C 

QRADSUM=0 0D00 
DO 30 KK=1,NRS 
QRAD (KK) = 0. 0D00 
DO 40 J = 1,NRS 

QRAD(KK) =QRAD(KK) + BBINV(KK,J)*DD(J) 

40 CONTINUE 

QRADSUM=QRADSUM+QRAD(KK) 

30 CONTINUE 
RETURN 
END 
C 

SUBROUTINE SHELL(C8,C9,C10,C11) 

CHARACTERS PHASE (40,20) ,NC 

REAL MDOT, KW, KS , KL, KLE, MT, MS, ML, K(40,20), KIP(40,20) , 
JcKIM(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU,NUL,KF 
COMMON DT,DRR,D2Z,DELW0,DELWI,DELWS,RI1,KIP,BBINV(18,18),QRAD(18), 
4tRI2,R01,R02,Q,TFI,MDOT,CPF,KW,CPW,RHOW,KS,CPS,RHOS,FF(18,18) l 
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&KL,CPL 1 RHOL,TM 1 HM I KLE,Z2,Z3,Z4,PI,VOLT,Q1 I Q2 I VOL(40,20),RCML(20), 

*MT,MS,VOLS,ML,VOLL,ESHSUM,EFSUM,RA,NU,CHARL,C7,C12 1 ECON(20) > 

&TV,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20),AREAOLD(3) I 

4tYF(40,20),DR(40),DZ(20),Z(20),R(40),EO(40,20),E(40,20),AREA(30) > 

&TFO (20) ,TF (20) ,C1 (40,20) ,C2(40,20) ,C3(40,20) ,C4(40,20),RCMLO(20), 
4tTFOUT,EE(40,20),TIMEPRl,TIMEPR2,TIMEP,TIMEF ,TFP .TFF.RV.RM, 
&TIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL,EEL(40,20) P 
&EELO(40,20),EEL02(40,20),PR,EML(40,20), 
tMFL.QF.II.IVl.lW.JJ.N.PHASE.NC, 

&ZKLE (1 0) ,ZN U (1 0) ,ZRA(1 0) ,ZCHARL(1 0) ,ZN UAVG 

UP-DATE CONDUCTORS AND DETERMINE ENERGY TRANSFER 


DO 100 J=2,JJ-1 

KIP(1 ,J) = K(2,J)*K(1 , J)*(DR(2) + DR(1))/(K(2,J)*DR(1) 
k +K(1,J)*DR(2)) 

KJP(1,J)=K(1,J + 1)*K(1,J)*(DZ(J+1)+DZ(J))/(K(1,J+1)*DZ(J) 

k +K(1,J)*DZ(J+1)) 

KJM(1 , J) = K(1,J— 1)*K(1 ,J)*(DZ(J-1)+DZ(J))/(K(1,J-1)*DZ(J) 

L +K(1,J)*DZ(J-1)) 

E(l,J)=EO(l,J)+Cl(l,J)*K!P(l,J)*(TO(2,J)-TO(l,J)) 
k +C8*(TFO(J)-TO(l,J)) 

k +C3(l,J)»KJP(l,J)*fTO(l,J+l)-TO(l,J)) 

k +C4(l,J)*KJM(l,J)*(TO(l,J— 1)— TO(l,J)) 

KIM(II,J)=K(II-1,J)*K(II,J)*(DR(II-1)+DR(II))/(K(II-1,J)*DR(II) 

k +K(II,J)*DR(II-1)) 

KJP(II,J) = K(II,J + 1)*K(II,J)*(DZ(J+1)+DZ(J))/(K(II,J+1)*DZ(J) 

k +K(II,J)*DZ(J + 1)) 

KJM(II,J) = K(II,J-1)*K(II,J)*(DZ(J-1)+DZ(J))/(K(I«,J-1)*DZ(J) 

k +K(II,J)*DZ(J-1)) 

E(II,J) = EO(II,J)+Q*C10-QRAD(J-1+JJ-2)*DT 
k +C2(ll,J)*KIM(ll,J)*(TO(ll-l,J)-TO(ll,J)) 

k +C3(ll,J)*KJP(ll,J)*(TO(ll,J+l)-TO(ll,J)) 

k +C4(ll,J)*KJM(ll,J)*(TO(ll,J-l)-TO(ll,J)) 

ECON(J)=C2(ll,J)*KIM(ll,J)*(TO(ll-l. J )- TO ( ll . J )) 

100 CONTINUE 
DO 105 1=2, IV1 

KIP(I,1)=K(I+1,1)*K(I,1)*(DR(I+1)+DR(I))/(K(I + 1,1)*DR(I) 

k +K(I,1)*DR(I + 1)) 

KIM(I,1)=K(I-1,1)*K(I,1)*(DR(I-1) + DR(I))/(K(I-1,1)*DR(I) 

k +K(I,1)*DR(I-1)) 

KJP(I,1)=K(I,2)*K(I,1)*(DZ(2)+DZ(1))/(K(I,2)*DZ(1) 
k +K(I,1)*DZ(2)) 

E(l,l) = EO(l,l)+Cl(l,l)*KIP(l,l)*(TO(l + l,l)-TO(l,l)) 
k +C2(l,l)*KIM(l,l)*(TO(l-l,l)-TO(l,l)) 

k +C3(l,l)*KJP(l,l)*(TO(l,2)-TO(l,l)) 

KIP(I,JJ)= K(l + 1 ,JJ)*K(I,JJ)*(DR(I + 1) + DR(I))/(K(I + 1,JJ)*DR(I) 
k +K(I,JJ)*DR(I+1)) 

KIM(I,JJ)=K(I-1,JJ)*K(I,JJ)*(DR(I-1) + DR(I))/(K(I-1,JJ)*DR(I) 
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* +K(I,JJ)*DR(I-1)) 

KJM(I,JJ)=K(I,JJ-1)*K(I,JJ)*(DZ(JJ-1)+DZ(JJ))/(K(I,JJ-1)»DZ(JJ) 

It +K(I,JJ)*DZ(JJ-1)) 

E(l,JJ)=EO(l,JJ)+Cl(l 1 JJ)*K!P(l I JJ)*(TO(H-l,JJ)-TO(l,JJ)) 

* +C2(I,JJ)*WM(I,JJ)*(T0(I-1 I JJ)-T0(I 1 JJ)) 

It +C4(I,JJ)*KJM(I,JJ)*(T0(I,JJ— 1)— TO(l,JJ)) 

105 CONTINUE 

JRANGE = 2*(JJ-2) + 1 
JJJ=1 

DO 106 l=IV2,ll-l 

WP(I I 1)=K(I+1 1 1)*K(I,1)*(DR(I+1)+DR(I))/(K(I+1 I 1)*DR(I) 

* +K(I,1)*DR(I+1)) 

WM(I,1)=K(I-1,1)*K(I,1)*(DR(I-1)+DR(I))/(K(I-1 ) 1)*DR(I) 

It +K(I,1)*DR(I-1)) 

KJP(I,1)=0 0D00 

E(l.l)=EO(l,l)+Cl(l,l)*KIP(l,l)*(TO(l+l,l)— TO(l,l)) 

It +C2(I,1)*KIM(I,1)*(TO(I-1 i 1)-TO(I i 1)) 

It +C3(l 1 l)*KJP(l,l)*frO(l 1 2)-TO(l,l)) 

It -QRAD(JRANGE)*DT*AREAOLD(JJJ)/AREA(JRANGE) 

KIP(I, JJ) = K(l + 1, JJ)*K(I, JJ)*(DR(I + 1) + DR(I))/(K(I+1,JJ)*DR(I) 

It +K(I,JJ)«DR(I + 1)) 

KIM(I,JJ)=K(I-1,JJ)*K(I,JJ)*(DR(I-1)+DR(I))/(K(I— 1,JJ)*DR(I) 

It +K(I,JJ)*DR(I-1)) 

KJM(l,JJ) = 0 ODOO 

E(l,JJ) = EO(l 1 JJ)+Cl(l > JJ)*KIP(l,JJ)*fTO(l+l,JJ)-TO(l 1 JJ)) 

It +C2(I,JJ)*KIM(I,JJ)*(T0(I-1,JJ)-T0(I,JJ)) 

It +C4(I 1 JJ)*KJM(I,JJ)*(T0(I,JJ-1)— TO(l,JJ)) 

It -QRAD(JRANGE+1)*DT*AREA0LD(JJJ)/AREA(JRANGE+1) 

JJJ=JJJ + 1 
106 CONTINUE 

KIP(1,1) = K(2 1 1)*K(1,1)*(DR(2) + DR(1))/(K(2,1)*DR(1) 

It +K(1,1)*DR(2)) 

KJP(1 , 1) = K(1 ,2) »K(1 , 1)»(DZ(2) + DZ(1))/(K(1 ,2) *DZ(1) 

It +K(1,1)*DZ(2)) 

E(1,1)=E0(1 I 1)+C1(1,1)*KIP(1 ) 1)*(T0(2,1)-T0(1 1 1)) 

It +C9*(TF0(1)-T0(1,1)) 

It +C3(l,l)*KJP(l l l)*(TO(l,2)-TO(l,l)) 

KIP(1,JJ) = K(2 1 JJ)*K(1 j JJ)*(DR(2)+DR(1))/(K(2,JJ)*DR(1) 

It +K(1,JJ)*DR(2)) 

KJM(1 , JJ) = K(1 , JJ — 1)*K(1 , JJ)*(DZ(JJ-1) + DZ(JJ))/(K(1 j JJ-1)*DZ(JJ) 

It +K(1,JJ)*DZ(JJ-1)) 

E(1 , JJ) = E0(1 , JJ) +C1 (1 , JJ)*KIP(1 , JJ)*(T0(2,JJ)-T0(1 ) JJ)) 

It -I- C9*(TF0(JJ) — T0(1,JJ)) 

It +C4(1,JJ)*KJM(1,JJ)*(T0(1 ) JJ— 1)— T0(1,JJ)) 

KIM(II 1 1)=K(II-1 1 1)*K(II,1)*(DR(II-1)+DR(II))/(K(II-1,1)»DR(II) 

It +K(II,1)*DR(II-1)) 

KJP(II i 1) = K(II i 2)*K(II,1)*(DZ(2) + DZ(1))/(K(II,2)*DZ(1) 

It +K(II,1)*DZ(2)) 

E(II,1)=E0(II,1)+Q*C11 
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* +C2(II,1)*KIM(II,1)*(T0{II-1,1)-T0(II,1)) 

Sc +C3(ll 1 l)*KJP(ll,l)*(TO(ll I 2)-TO(ll,l)) 

KIM(II,JJ) = K(II-1 i JJ)*K(II j JJ)*(DR(H-1)+DR(II))/(K(II-1 i JJ)* 

Sc DR(II) + K(II,JJ)*DR(II-1)) 

KJM(II,JJ)=K(II,JJ-1)*K(II,JJ)*(DZ(JJ-1) + DZ(JJ))/(K(II,JJ— 1)* 

Sc DZ(JJ)+K(ll,JJ)*DZ(JJ-lj) 

E(ll,JJ) = EO(ll l JJ)+Q*Cll 

Sc +C 2 (ll,JJ)*WM(ll,JJ)*fTO(ll-l,JJ)-TO(ll,JJ)) 

Sc +C4(ll ) JJ)*KJM(ll,JJ)*(TO(ll 1 JJ— 1)— TO(ll,JJ)) 

DETERMINE ELEMENT TEMPERATURES 

DO 110 J=1,JJ 

EE(l,J) = E(l,J)/(RHOW*VOL(l,J)) 

T(1,J) = (EE(1,J)/CPW)+TM 
EE(ll,J) = E(ll l J)/(RHOW*VOL(ll,J)) 

T(II,J) = (EE(II,J)/CPW) +TM 
110 CONTINUE 
DO 115 1=2,11-1 
EE(l,l)=E(l,l)/(RHOW*VOL(l,l)) 

T(l,l) = (EE(I,1)/CPW) +TM 
EE(l,JJ)=E(l 1 JJ)/(RHOW*VOL(l,JJ)) 

T(I,JJ) = (EE(I,JJ)/CPW)+TM 
115 CONTINUE 
RETURN 
END 
C 

SUBROUTINE FLUID(C13,C14,C15,C16,C17) 

CHARACTERM PHASE(40,20),NC 

REAL MDOT,KW,KS,KL,KLE,MT,MS,ML 1 K(40 ) 20),K1P(40 1 20) ) 
&KIM(40,20),KJP(40,20),KJM(40 1 20) 1 MFL,MUF,NU 1 NUL 1 KF 

COMMON DT.DRR.DZZ.DELWO.DELWI.DELWS.Rll.KIP.BBINVaMfyQRADflS), 
fcRtt.ROl.ROZ.Q.TFI.MDOT.CPF.KW.CPW.RHOV^KS.CPS.RHOS.FFllMS), 
AKL,CPL,RHOL,TM,HM 1 KLE,Z2,Z3,Z4,PI,VOLT,Q1,Q2 I VOL(40 J 20) 1 RCML(20) 1 
&MT,MS I VOLS,ML,VOLL I ESHSUM J EFSUM,RA,NU,CHARL 1 C7 1 C12 1 ECON(20) 1 
4£TV,TO(40 1 20),T(40 I 20) ) K,RHO(40,20) 1 XF(40 I 20) I TFIN(20),AREAOLD(3), 
4sYF(40,20),DR(40) 1 DZ(20) ) Z(20),R(40) 1 EO(40,20),E(40,20) ) AREA(30) l 
&TFO(20),TF(20),C1(40,20) ) C2(40,20) 1 C3(40 1 20) 1 C4(40 1 20),RCMLO(20), 
&TFOUT,EE(40 1 20) 1 TIMEPR1,TIMEPR2 1 T1MEP 1 TIMEF,TFP,TFF 1 RV 1 RM, 
&TIME2,TIME3 ) TIME4,TF2,TF3 1 TF4,EINIT 1 ECAN,EBAL,EEL(40 ) 20) 1 
4tEELO(40,20),EEL02(40,20),PR 1 EML(40 ) 20) 1 
AMFL.QF.II.IVIJW^J.N.PHASE.NC, 

1ZKLE (1 0) ,ZN U (10) ,ZRA(10) ,ZCHARL(1 0) ,ZN UAVG 
C 

C DETERMINE FLUID TEMPERATURE DISTRIBUTION BASED ON 
C A "PSEUDO-STEADY STATE" ASSUMPTION 
C 

TF(1) = (C13*T(1,1)+C14*TFI)/C15 
TFI N (2) = 2 ODOO*TF (1 ) -TF I 
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DO 200 J=2,JJ— 1 

TF(J) = (C16*T(1,J)+C14*TFIN(J))/C17 
TFIN(J+1)=2 0D00*TF(J)-TFIN(J) 

200 CONTINUE 

TF(JJ)=(C13*T(1 1 JJ)+C14*TFIN(JJ))/C15 
TFOUT =2.0D00*TF(JJ)-TFIN(JJ) 

QF=MDOT*CPF*(TFOUT-TFI) 

EFSUM = EFSUM + QF*DT 
RETURN 
END 
C 

SUBROUTINE SALT(DUM) 

CHARACTERM PHASE(40,20),NC 
REAL MDOT 1 KW 1 KS 1 KL,KLE,MT,MS J ML,K(40 1 20) 1 KIP(40 I 20) ) 
4tKIM(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU,NUL,KF,EEO(40,20) 

COMMON DT,DRR,DZZ,DELW0,DELWI,DELWS,RI1,KIP,BBINV(18,18),QRAD(18), 
leRI2,R01,R02 1 Q,TFI,MD0T,CPF,KW,CPW,RH0W,KS,CPS,RH0S,FF(18,18), 
tKL,CPL 1 RHOL,TM 1 HM,KLE 1 Z2,Z3,Z4,Pl 1 VOLT p Ql 1 Q2 I VOL(40,20) l RCML(20) 1 
4tMT,MS,VOLS, ML, VOLL,ESHSUM, EFSUM, RA,NU,CHARL,C7,C12,ECON(20) ( 
4cTV,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20) 1 AREAOLD(3), 

8iYF (40,20), DR(40),DZ(20),Z(20),R(40),EO(40, 20), E(40,20),AREA(30), 

*TFO(20),TF(20),C1(40 I 20),C2(40,20),C3(40,20),C4(40,20),RCMLO(20), 

*TFOUT,EE(40,20),TIMEPR1,TIMEPR2,TIMEP,TIMEF,TFP,TFF,RV,RM, 

*TIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL,EEL(40,20), 

*EELO(40,20) I EEL02(40,20),PR,EML(40,20), 

AMFL,QF, II, IV1,IV2,JJ,N, PHASE, NC, 
*ZKLE(10),ZNU(10),ZRA(10),ZCHARL(10),ZNUAVG 
C 

C UP-DATE CONDUCTORS, DETERMINE ENERGY TRANSFER, AND 
C UP-DATE PCM ELEMENT TEMPERATURES AND PROPERTIES 
C 

ML=00D00 
DO 300 I=2,IV1-1 
DO 305 J =2, JJ — 1 

KIP(I, J) = K(l + 1 ,J)*K(I,J)*(DR(I + 1) + DR(I))/(K(I+ 1, J)*DR(I) 

L +K(I,J)*DR(I + 1)) 

KIM(I 1 J)=K(I-1,J)*K(I,J)*(DR(I-1)+DR(I))/(K(I-1,J)*DR(I) 
k + K(I,J)*DR(I-1)) 

KJP(I,J)=K(I,J+1)*K(I,J)*(DZ(J+1)+DZ(J))/(K(I,J+1)*DZ(J) 
k +K(I,J)*DZ(J+1)) 

KJM(I,J) = K(I,J-1)*K(I,J)*(DZ(J-1)+DZ(J))/(K(I,J-1)*DZ(J) 
k +K(I,J)*DZ(J-1)) 

E(l,J)=EO(l,J)+Cl(l,J)*KIP(l,J)*(TO(l+l,J)— TO(l,J)) 
k +C2(I,J)*KIM(I,J)*(TO(I-1 i J)-TO(I,J)) 

* +C3(l,J)*KJP(l,J)*(TO(l,J+l)-TO(l,J)) 

k +C4(l,J)*KJM(l,J)*(TO(l,J-l)-TO(l,J)) 

EE(l,J)=E(l,J)/(RHO(l,J)»VOL(l,J)) 

IF (EE(I,J) GT HM) GO TO 310 
IF (EE(l,J).LT.O D00) GO TO 315 
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T(I,J)=TM 

XF(I,J)=EE(I,J)/HM 

YF(I,J) = 1 DOO/(1.DOO+(RHOL/RHOS)*(1.DOO/XF(I,J)— l.DOO)) 

K(l, J) = (1 DOO-XF(l,J))*KS+XF(l, J)*ZKLE(J) 

RHO(I,J) = (1.DOO-YF(I,J))*RHOS+YF(I,J)*RHOL 
PHASE(I, J) = ’MUSH’ 

GO TO 304 

310 T(I,J) = (EE(I,J)-HM)/CPL+TM 
K(I,J)=ZKLE(J) 

RHO(l,J)=RHOL 
XF(I, J) = 1.0D00 
PHASE(I,J) = ’UQ’ 

GO TO 304 

315 T(I,J) = (EE(I,J)/CPS)+TM 
K(I,J) = KS 
RHO(l,J)=RHOS 
XF(I,J)=O.ODOO 
PHASE(l,J) = ’SOL’ 

304 ML=ML+XF(l,J)*RHOL*VOL(l,J) 

305 CONTINUE 
300 CONTINUE 

I=IV1 

DO 320 J=2,JJ-1 

WM(I j J) = K(I-1 1 J)*K(I,J)*(DR(I-1)+DR(I))/(K(I-1,J)*DR(I) 

L +K(I,J)*DR(I-1)) 

KJP(I,J)=K(I i J+1)*K(I,J)*(DZ(J + 1)+DZ(J))/(K(I 1 J+1)*DZ(J) 

Sc +K(I,J)*DZ(J + 1)) 

KJM(I,J) = K(I,J-1)*K(I,J)*(DZ(J-1) + DZ(J))/(K(I 1 J— 1)*DZ(J) 

Sc +K(I,J)*DZ(J-1)) 

E(l,J) = EO(l,J)-ECON(J) 

Sc +C2(l,J)*KIM(l 1 J)*(TO(l-l 1 J)-TO(l,J)) 

Sc +C3(l 1 J)*KJP(l,J)*(TO(l,J + l)-TO(l,J)) 

Sc -F C4(l, J)*KJM(l,J)*(TO(l, J — 1) — TO(l, J)) 

Sc -QRAD(J-1)*DT 

EE(I, J) = E(l,J)/(RHO(l,J)*VOL(l, J)) 

IF (EE(I,J).GT HM) GO TO 330 
IF (EE(I,J) LT.O D00) GO TO 335 
T(I,J)=TM 
XF(I,J)=EE(I 1 J)/HM 

YF (I , J) = 1 D00/(1 . D00 + (RHOL/RHOS) * (1 DOO/XF (I , J) —1 . D00)) 
K(I,J) = (1 D00-XF(l,J))*KS+XF(l, J)*ZKLE(J) 

RHO(I,J) = (1.DOO-YF(I 1 J))*RHOS+YF(I,J)*RHOL 
PHASE(I,J) = ’MUSH’ 

GO TO 324 

330 T(I,J) = (EE(I,J)-HM)/CPL+TM 

k(i,J)=zkle(j) 

RHO(l,J)=RHOL 
XF(l,J) = l 0D00 
PHASE(I,J) = ’UQ’ 
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GO TO 324 

335 T(l, J) = (EE(I,J)/CPS) +TM 
K(I,J)=KS 
RHO(l,J)=RHOS 
XF(l,J)=0 0D00 
PHASE(l,J)= , SOL’ 

324 ML=ML+XF(l,J)*RHOL*VOL(l,J) 

320 CONTINUE 
MFL=ML/MT 
RETURN 
END 
C 

SUBROUTINE UPDATE(N1,N2,N3,WF) 

CHARACTERM PHASE(40,20),NC 

REAL MDOT,KW,KS,KL,KLE,MT I MS,ML,K(40 > 20),KIP(40 I 20), 
&KIM(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU,NUL,KF 
COMMON DT,DRR,DZZ,DELW0,DELWI,DELWS,RI1,KIP,BBINV(18,18),QRAD(18), 
tRI2 1 R01,R02 1 Q 1 TFI 1 MD0T,CPF > KW,CPW I RH0W l KS l CPS,RH0S,FF(18,18) 1 
4tKL J CPL ) RHOL 1 TM,HM l KLE,Z2 1 Z3,Z4,Pl 1 VOLT,Ql J Q2 1 VOL(40 1 20),RCML(20) 1 
*MT,MS,VOLS,ML,VOLL,ESHSUM,EFSUM,RA,NU,CHARL,C7,C12,ECON(20), 
4tTV,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20) 1 AREAOLD(3), 
&YF(40,20),DR(40),DZ(20),Z(20),R(40),EO(40,20),E(40,20),AREA(30), 

4eTFO(20),TF (20), Cl (40,201,02(40,20), C3(40,20),C4(40,20) J RCMLO(20) ) 
&TFOUT,EE(40,20),TIMEPR1,TIMEPR2,TIMEP,TIMEF I TFP,TFF,RV,RM, 
4tTIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL,EEL(40 l 20), 
&EELO(40,20),EEL02(40,20),PR,EML(40,20), 

4cMFL,QF, II, IV1,IV2,JJ,N, PHASE, NC, 

1ZKLE (1 0) ,ZN U (1 0) ,ZRA(1 0) ,ZCH ARL(1 0) , ZN UAVG 
TIMEM=(N-l)*DT/60 0D00 
IDIFFO= (II — 1) — IV1 
IDIFFN = (II — 1 j — IV1 

VOID VIEW FACTORS 

IF (IDIFFN NE IDIFFO) CALL VOIDFF(ES,EW) 

SALT NATURAL CONVECTION CORRELATION 
IF(NC EQ ’ON ’) CALL CONV(DUM) 

GLOBAL HEAT TRANSFER AND ENERGY BALANCE 
CALL ENERGY(EBALP,TIMEM) 

DETERMINE SIDE WALL HEAT TRANSFER FRACTIONS 
CALL WALLF RAC (WALLF R1 ,WALLFR2,WALLFR3) 

PRINT RESULTS TO OUTPUT FILES 
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C 

CALL 0 UTP UT(TI M E M , N 1 , N2, WALLF R1 , WALLFR2 , WALLFR3.EBALP, WF) 

UP-DATE INPUT HEAT FLUX 

IFfTIMEM.GE 54 63) Q=Q2 
IF(TIMEM GE 72.8) Q=Q3 

UP-DATE INPUT COOLING FLUID INLET TEMPERATURE 

IFfTIMEM GE TIME2.AND.TIMEM.LT.TIME3) GO TO 920 
IFfTIMEM GE TIME3) GO TO 930 
GO TO 905 
920 TIMEP=TIME2 
TIMEF=TIME3 
TFP=TF2 
TFF=TF3 
GO TO 905 
930 TIMEP=TIME3 
TIMEF=TIME4 
TFP=TF3 
TFF=TF4 

905 TFI=TFP+(TIMEM-TIMEP)*(fTFF-TFP)/(TIMEF— TIMEP)) 

950 RETURN 
END 
C 

SUBROUTINE VOIDFF(ES.EW) 

CHARACTERS PHASE(40,20),NC 

REALM DOT, KW, KS , KL, KLE , MT, MS , ML, K(40,20) , KIP (40,20) , 
tKIM(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU,NUL,KF, 
tBB(18,18),EPS(30),FFSUM(18),F12(18),F13(18), 

&F23(18),F34(18) 

COMMON DT.DRR.DZZ.DELWO.DELWI.DELWS.Rll.KIP.BBINVflMSJ.QRADfiS), 

tRI2.R01.R02, Q.TFI.MDOT, CPF, KW.CPW.RHOW.KS.CPS.RHOS.FFflS.lS), 
tKL,CPL,RHOL,TM,HM,KLE,Z2,Z3,Z4,PI,VOLT,Ql,Q2,VOL(40,20),RCML(20), 
tMT,MS 1 VOLS l ML,VOLL l ESHSUM,EFSUM,RA,NU,CHARL,C7,C12,ECON(20), 
tTV,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20),AREAOLD(3), 
&YF(40 1 20),DR(40),DZ(20),Z(20),R(40),EO(40,20),E(40,20),AREAp0), 
tTFO(20),TF(20),Cl(40,20),C2(40,20),C3(40,20),C4(40,20),RCMLO(20), 
tTFOUT,EE(40, 20), TIMEPR1.TIMEPR2, TIMEP, TIMEF.TFP.TFF.RV.RM, 
tTIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL,EEL(40,20), 
tEELO(40,20),EEL02(40,20),PR,EML(40,20), 
tMFL.QF, II, IV1,IV2,JJ,N, PHASE, NC, 
tZKLE(10),ZNU(10),ZRA(10),ZCHARL(10),ZNUAVG 
NRS =2*(J J — 2) +2*(II-1-IV1) 

DO 10 J — 1 , JJ — 2 
EPS(J) = ES 
10 CONTINUE 

DO 20 J=JJ-1,NRS 
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EPS(J)=EW 
20 CONTINUE 
DO 30 J=1,JJ— 2 
AREA(J) =2 0D00*PI*RV*DZZ 
AREA(J + JJ-2) =2 0D00*PI*ROl*DZZ 
30 CONTINUE 

JRANGE=2*(JJ-2) + l 
l = IV2 
JJJ = 1 

DO 40 J=JRANGE,NRS-1,2 
AREA(J) =2 ODOO*PI*R(I)*DR(I) 

AREAOLD(JJJ) =AREA(J) 

AREA(J + 1) = AREA(J) 

1 = 1+1 

JJJ=JJJ+1 

40 CONTINUE 

DO 11 KK=2*(JJ-2)+3,NRS-l,2 
AREA(JRANGE) = AREA(JRANGE) + AREA(KK) 

AREA(JRANGE + 1) = AREA(JRANGE) 

11 CONTINUE 
NRS=2*(JJ-2)+2 
LDBBINV=NRS 
LDBB=NRS 

DEFINE THE FF(KK,J) VIEW FACTOR MATRIX FOR SURFACE 
C ELEMENTS COMPRISING THE ANNULAR GEOMETRY VOID WITH 
C AN INNER SURFACE "1", AN OUTER SURFACE "2", AND 
C SIDE SURFACES "3" AND "4" 

C 

DZZ1=DZZ 
DO 41 KK=l,JJ-2 
DO 42 J=JJ-l,2*(JJ-2) 

ZKJ=DZZl*(ABS(J-(JJ-2)-KK)-l) 

IF (ZKJ LT 0 0D00) GO TO 42 

FF(KK,J) = FF1K2J(RV,R01 1 DZZ1 1 DZZ1,ZKJ,PI) 

FF(J,KK)=FF(KK,J)*AREA(KK)/AREA(J) 

42 CONTINUE 

FF (KK+ J J-2, KK) = FF21 (RV, ROl , DZZ1 ,PI) 

FF(KK,KK+ JJ-2) =FF(KK+ JJ-2, KK)*AREA(KK+JJ-2)/AREA(KK) 

41 CONTINUE 

DO 43 KK=l,JJ-2 
l=IV2 

ZKJ=DZZ1*KK 

DO 44 J=2*(JJ-2) + l,NRS-l,2 
DR2=R01-RV 
R2=(R01 + RV)/2 D00 
FF(KK I J) = FF1K3J(RV,R2,DZZ1,DR2,ZKJ,PI) 

FF(J,KK) = FF(KK,J)*AREA(KK)/AREA(J) 

FF(JJ— 1-KK,J+1)=FF(KK,J) 
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FF(J+1,JJ-1-KK) = FF(J,KK) 

1=1 + 1 

44 CONTINUE 
43 CONTINUE 

DZZ1=DZZ 

FF(JJ-1,JJ-1)=FF22(RV I R01 I DZZ1,PI) 

FFJJ1=FF(JJ,1) 

DO 45 KK= JJ — 1 ,2*(JJ — 2) 

DO 46 J= JJ — 1,2* (JJ— 2) 

IF(KK.EQ.J) GO TO 46 

FF(KK,J)=FF2K2J(RV I R01,DZZ1,KK,J,PI,FFJJ1) 

46 CONTINUE 
FF(KK,KK) = FF(JJ-1 I JJ-1) 

45 CONTINUE 

DO 47 KK=JJ-l,2*(JJ-2) 
l=IV2 

DO 48 J=2*(JJ-2) + l,NRS-l,2 
ZKJ = DZZ1*(KK — (JJ — 2)) 

R1 = (RV+R01)/2 D00 
DR1=R01-RV 

FF(KK,J)=FF2K3J(R1,R01,DZZ1 J DR1,ZKJ,PI) 

FF(J,KK)=FF(KK,J)*AREA(KK)/AREA(J) 

FF(3*(JJ — 2) + 1-KK, J+ 1) = FF(KK,J) 

FF(J+l,3*(JJ-2) + l-KK)=FF(J,KK) 

1 = 1 + 1 

48 CONTINUE 

47 CONTINUE 
ZKJ=Z3-Z2 
13= IV2 
1=13 

JJJ=2*(JJ — 2) +2 
KKK=2*(JJ-2) + l 
DO 49 KK=KKK,NRS-1,2 
Rl = (ROl + RV)/2.D00 
DR1 = R01-RV 
DO 50 J=JJJ,NRS,2 
DR2=R01-RV 
R2= (ROl + RV)/2 D00 
FF(KK I J)=FF3K4J(R1,R2 1 DR1,DR2 I ZKJ J PI) 
FF(J,KK)=FF(KK,J)*AREA(KK)/AREA(J) 

1=1 + 1 

50 CONTINUE 
13=13+1 
1 = 13 

JJJ = JJJ +2 
KKK=KKK+2 

49 CONTINUE 
C 

C DEFINE "BB" MATRIX SUCH THAT {BB}{QRAD}={DD> WHERE {BB> IS 


non 
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C THE 'VIEW FACTOR-EMITTANCE" MATRIX AND {DD> IS THE 
C TEMPERATURE-TO-THE-FOURTH-POWER DIFFERENCE" MATRIX 
C 

NRS=2*(JJ-2)+2 
DO 55 KK=1,NRS 
DO 60 J=1,NRS 
DELTA=0 0D00 
IF(KKEQ J) DELTA=1.0D00 

BB(KK,J)=(DELTA/EPS(J)-FF(KK,J)*(1.0D00— EPS(J))/EPS(J))/AREA(J) 
60 CONTINUE 
55 CONTINUE 

INVERT {BB> MATRIX USING IMSL ROUTINE DUNRG 

CALL DLINRG(NRS,BB,LDBB,BBINV,LDBBINV) 

QRADSUM = 0 0D00 
DO 70 KK=1,NRS 
FFSUM (KK) = 0 ODOO 
DO 80 J = 1,NRS 

FFSUM(KK)=FFSUM(KK) + FF(KK,J) 

80 CONTINUE 

WRITE (14,1001) KK,FFSUM(KK),QRAD(KK) 

QRADSUM =QRADSUM+QRAD(KK) 

70 CONTINUE 

WRITE (14,1004) QRADSUM 

WRITE (14,1002) ((FF(KK,J),J=1,NRS),KK=1,NRS) 

1001 F0RMAT(I4,2X,F8 6,2X,F9 5) 

1002 FORMAT(18(F6 4, IX)) 

1004 FORMAT(16X,F9.5) 

RETURN 

END 

C 

FUNCTION FF13(R113,R213,L13,PI) 

REAL L13 
RR=R213/R113 
A=L13**2+R213**2-R1 13**2 
B = L13**2-R213**2 + R1 13**2 
C = L13**2+ R213**2+ R1 13**2 

FF13=(1 DOO/(2.DOO*PI))*(DARCOS(B/A)-(R113/(2.DOO*L13))* 
tt (DSQRT(C**2/R113**4-4D00*RR**2)*DARCOS(B/(RR*A))+ 

* (B/R113**2)*DARSIN(l.D00/RR)-(PI/2 D00)*(A/R113**2))) 

RETURN 

END 

C 

FUNCTION FF21(R121,R221,L21F,PI) 

REAL L21.L21F 
RR=R221/R121 
L21=L21F/R121 
A=L21**2+RR**2-1 D00 
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B=L21**2-RR**2+1 DOO 

FF21 = l.DOO/RR-(l DOO/(PI*RR))*(DARCOS(B/A)-(1.DOO/(2.DOO*L21))* 

& (DSQRT ((A+ 2. DOO)**2-4.DOO*RR**2)*DARCOS(B/(RR*A)) + 

St B*DARSIN(1 . DOO/RR) -PI*A/2. DOO)) 

RETURN 

END 

FUNCTION FF22(R122,R222,L22F,PI) 

REAL L22.L22F 
RR-R222/R122 
L22=L22F/R122 
A=L22**2+RR**2— l.DOO 
B = L22**2— RR**2 + 1 . DOO 

FF22=1.D00-1.D00/RR+(2.D00/(PI*RR))*DATAN(2.D00*DSQRT(RR**2— l.DOO) 
L /L22)-(L22/(2DOO*PI*RR))*((DSQRT(4DOO*RR**2+L22**2)/L22)* 

L DARSIN((4 DOO*(RR**2-l.DOO) + (L22**2/RR**2)*(RR**2— 2.DOO))/ 

L (L22**2+ 4 DOO* (RR**2-1 DOO))) -DARSIN ((RR**2-2.DOO)/RR**2) + 

L (PI/2 DOO) * (DSQRT (4. D00*RR**2 + L22**2)/L22-l DOO)) 

RETURN 

END 

: 

FUNCTION FF34(R134 J R234,L34,PI) 

REAL L34 

FF34=1.D00— (2 D00*R134*L34/(R234**2— R134**2))* 

St FF13(R134,R234 1 L34,PI)- 

St (2 D00*R234*L34/(R234**2— R134**2))*(1.D00— 

St FF22(R134,R234,L34,PI)— 

St FF21(R134,R234,L34,PI))/2.DOO 
RETURN 
END 
C 

FUNCTION FF1 K2J(R1 ,R2,L > Y,D,PI) 

REAL L 

IF (D EQ O ODOO) FF1K2J = ((L+D)/L)*FF13(R1,R2,L+D,PI) + 

St ((Y+D)/L)*FF13(R1 I R2,Y+D,PI)-((L+D+Y)/L)*FF13(R1,R2 I L+D+Y,PI) 

IF (D EQ 0 ODOO) GO TO 12 

FF1K2J = ((L+D)/L)*FF13(R1 1 R2,L+D 1 PI) + ((Y+D)/L)*FF13(R1,R2,Y+D,PI)- 
St (D/L)*FF13(R1,R2,D,PI)-((L+D+Y)/L)*FF13(R1,R2,L+D+Y,PI) 

12 RETURN 
END 
C 

FUNCTION FFUG-HRI^DZICDRJ.ZKJ.PI) 

REAL L 

R2P=R2+DRJ/2DOO 
R2M=R2-DRJ/2.DOO 
RDIFF = DABS(R2-R1) 

ZDIFF=DABS(ZKJ-DZK) 

IF (RDIFF.GT.DRJ) GO TO 39 
IF (ZDIFF GT.O.OIDOO) GO TO 23 
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FF1K3J=FF13(R1,R2P,ZKJ,PI) 

GO TO 19 

23 FF1K3J=(ZKJ/DZK)*FF13(R1,R2P I ZKJ,PI)- 

A ((ZKJ-DZK)/DZK)*FF13(R1,R2P,ZKJ-DZK,PI) 

GO TO 19 

39 IF (ZDIFF.GT.O 01D00) GO TO 24 

FF1K3J=FF13(R1 I R2P 1 ZKJ,PI)-FF13(R1 I R2M,ZKJ,PI) 

GO TO 19 

24 FF1K3J=(ZKJ/DZK)*(FF13(R1 1 R2P,ZKJ I PI)-FF13(R1, 

A R2M,ZKJ 1 PI))-((ZKJ-DZK)/DZK)*(FF13(R1, 

A R2P 1 ZKJ-DZK,PI)-FF13(R1,R2M,ZKJ-DZK,PI)) 

19 RETURN 
END 
C 

FUNCTION FF2K2J(R1,R2,ZZZ22,KK,J,PI,FFJJ1) 

REAL L 

ZERO=OODOO 

DIFFKJ = ABS(KK-J) + 1 D00 

DIFFKJ1=ABS(KK-J)*1.D00 

DIFFKJ2=DIFFKJ1-1 D00 

DIFFKJ3=DIFFKJ2-1.D00 

IF (DIFFKJ1 EQ 1 D00) GO TO 15 

IF (DIFFKJ1 EQ 2.D00) GO TO 16 

FF2K2J=DIFFKJ1*(1 D00-FF22(R1,R2,DIFFKJ1*ZZZ22 1 PI)— FF21(R1, 
A R2,DIFFKJ1*ZZZ22,PI))- 

St DIFFKJ2*(1 D00-FF22(R1,R2,DIFFKJ2*ZZZ22,PI)-FF21(R1, 

L R2 1 DIFFKJ2*ZZZ22,PI))/2 D00- 

* DIFFKJ*(1 D00— FF22(R1,R2,DIFFKJ*ZZZ22,PI)— FF21(R1, 

A R2,DIFFKJ*ZZZ22,PI))/2 D00- 

A FF1K2J(R1 1 R2,DIFFKJ1*ZZZ22,ZZZ22,ZER0 P PI)* 

A DIFFKJ1*R1/R2+ 

A FF1K2J(R1 I R2,DIFFKJ2*ZZZ22,ZZZ22 1 ZER0,PI)* 

A DIFFKJ2*R1/R2 

GO TO 17 

15 FF2K2J=(1 D00-FF22(R1,R2 J ZZZ22,PI)— FF21(R1,R2 1 ZZZ22,PI))— 

A (1 D00-FF22(R1,R2,2 D00*ZZZ22 1 PI)-FF21(R1, 

A R2,2 D00*ZZZ22,PI))-FFJJ1 

GO TO 17 

16 FF2K2J=DIFFKJ1*(1 D00-FF22(R1,R2,DIFFKJ1*ZZZ22 1 PI)— FF21(R1, 

A R2,DIFFKJ1*ZZZ22,PI))- 

A DIFFKJ2*(1 D00-FF22(R1 , R2, DIFFKJ2*ZZZ22, PI)— FF21 (R1 , 

A R2,DIFFKJ2‘ZZZ22,PI))/2 D00- 

A DIFFKJ*(1.D00-FF22(R1,R2,DIFFKJ*ZZZ22,PI)— FF21(R1, 

A R2,DIFFKJ*ZZZ22,PI))/2.DOO- 

A FF1K2J(R1,R2,DIFFKJ1*ZZZ22,ZZZ22,DIFFKJ3*ZZZ22,PI)* 

A DIFFKJ1*R1/R2+ 

A FF1K2J(R1,R2,DIFFKJ2*ZZZ22,ZZZ22,DIFFKJ3*ZZZ22,PI)* 

A DIFFKJ2*R1/R2 

17 RETURN 
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C 

FUNCTION FF2K3J(R1,R2,ZZZ 1 DR1,ZKJ,PI) 

REAL L 

RlP = Rl + DRl/2 DOO 
RlM=Rl-DRl/2 DOO 
RDIFF =DABS(R2-R1) 

DIFFKJ=ZKJ/ZZZ 
DIFFKJl = DIFFKJ-l.DOO 
IF (DIFFKJ EQ 1 DOO) GO TO 25 
IF (RDIFF. LT.DR1) GO TO 21 

FF2K3J = DIFFKJ*(1 D00-FF22(R1M,R2,DIFFKJ*ZZZ,PI)— 
k FF21(R1M,R2,DIFFKJ*ZZZ,PI))/2D00- 

k DIFFKJ1*(1.D00-FF22(R1M,R2 1 DIFFKJ1*ZZZ,PI)— 

k FF21(R1M,R2,DIFFKJ1*ZZZ 1 PI))/2.D00- 

k DIFFKJ*(1.D00-FF22(R1P,R2,DIFFKJ*ZZZ,PI)- 

k FF21(R1P,R2,DIFFKJ*ZZZ,PI))/2.D00+ 

k DIFFKJ1*(1.D00— FF22(R1P,R2 I DIFFKJ1*ZZZ,PI)— 

k FF21(R1P,R2,DIFFKJ1*ZZZ,PI))/2.D00 

GO TO 27 

21 FF2K3J=DIFFKJ*(1 D00-FF22(R1M,R2,DIFFKJ*ZZZ,PI)— 

k FF21(R1M,R2,DIFFKJ*ZZZ,PI))/2.D00- 

k DIFFKJ1*(1 D00-FF22(R1M,R2,DIFFKJ1*ZZZ 1 PI)- 

k FF21(R1M,R2,DIFFKJ1*ZZZ,PI))/2.D00 

GO TO 27 

25 IF (RDIFF LT DR1) GO TO 22 

FF2K3J = (1 D00-FF22(R1M,R2 1 ZZZ,PI)- 
k FF21(R1M,R2,ZZZ,PI))/2.D00- 

k (1 D00-FF22(R1P,R2,ZZZ,PI)- 

k FF21 (R1P 1 R2,ZZZ,PI))/2.D00 

GO TO 27 

22 FF2K3J = (1 D00-FF22(R1M,R2,ZZZ,PI)- 

k FF21(R1M,R2 ) ZZZ,PI))/2.D00 

27 RETURN 
END 
C 

FUNCTION FF3K4J(R1 1 R2,DR1,DR2,L 1 PI) 

REAL L 

R J P = R2 + D R2/2 DOO 
RJM = R2-DR2/2 DOO 
RKP = R1 +DR1/2 DOO 
RKM = Rl-DRl/2 DOO 
DIFFR1 = DABS(R2-R1) 

DIFFR2=DABS(RJM-RKP) 

IF (DIFFR1 GT 1 0D-05) GO TO 11 
FF3K4J = FF34(RKM,RJP,L,PI) 

GO TO 337 

11 IF (DIFFR2 GT.1.0D-05) GO TO 12 

FF3K4 J = (0 5D00) * (AJ DAK*F F34(RKM ,RJP, L, PI) -AAAK* 
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* FF34(RJM,RJP,L,PI)-FF34(RKM,RKP,L,PI)) 

FDA=FF3K4J 
GO TO 337 

12 FF3K4J= (0 5D00)*(AJCDAK*FF34(RKM,RJP,L,PI)-AABAK* 
k FF34(RKP,RJP,L,PI)-2.D00*FDA-FF34(RKM,RKP,L,PI)) 

337 RETURN 
END 
C 

SUBROUTINE CONV(DUM) 

CHARACTER‘4 PHASE(40,20),NC 

REAL MDOT, KW, KS, KL, KLE, MT,MS, ML, K(40,20) , WP(40,20), 
8tKlM(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU,NUL,KF,ZML(10) 

COMMON DT,DRR,DZZ,DELW0,DELWI,DELWS,RI1,KIP,BBINV(18,18),QRAD(18), 
&RI2 1 R01 > R02 I Q,TFI 1 MDOT 1 CPF,KW,CPW,RHOW,KS I CPS,RHOS I FF(18 I 18) I 
8tKL l CPL 1 RHOL,TM,HM 1 KLE 1 Z2,Z3,Z4,PI,VOLT ) Ql,Q2,VOL(40 1 20) l RCML(20), 
&MT,MS,VOLS,ML,VOLL,ESHSUM,EFSUM,RA,NU,CHARL,C7,C12,ECON(20), 
&TV,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20),AREAOLD(3), 
&YF(40,20),DR(40),DZ(20),Z(20),R(40),EO(40,20),E(40,20),AREA(30) I 
*TFO(20),TF(20),C1(40,20),C2(40,20),C3(40,20),C4(40,20),RCMLO(20), 
<tTFOUT,EE(40,20) l TIMEPRl ) TIMEPR2,TIMEP,TIMEF 1 TFP,TFF ) RV,RM, 
&TIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL,EEL(40,20) > 
&EELO(40,20),EEL02(40,20),PR,EML(40,20), 
fcMFL.QF.II.IVl.M.JJ.N, PHASE, NC, 
&ZKLE(10),ZNU(10),ZRA(10),ZCHARL(10),ZNUAVG 
MS=MT-ML 

RM = DSQRT(1 ODOO*RI2*RI2+MS/(RHOS*PI*(Z3-Z2))) 

CHARL=RV-RM 

IF (CHARL LE 0 0D00) GO TO 530 

C18= (Z3-Z2)/CHARL 

I=IV1 

IF (RV EQ.ROl) 1 = 11 
TV=0 0D00 
DO 500 J=2,JJ-1 
TV=TV+T(I,J) 

500 CONTINUE 
TV=TV/(JJ-2) 

C 

C NOTE FOR CASES WITHOUT VOID, I — > II AND TV=TSHELL C O.D. 

C 

RA=C7*(TV-TM)*CHARL**3 
IF (RA LE 0 0D00) GO TO 530 

NU=0 42D00*PR**0 O12DO0*RA**0 25D00*C18"(-° 3D00) 

IF (C18 LT 10 0D00) NU=0 22DOO*(RA*PR/(0.2DOO+PR))**0.28DOO 
k *C18**(-0.25D00) 

IF (NU LT 1 0D00) GO TO 530 
GO TO 531 

530 NU = 1 0D00 

531 KLE=KL*NU 
C 
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FREE CONVECTION WITH AXIAL DEPENDENCE 


ZMT=MT/(JJ-2) 

DO 540 J=2,JJ-1 
ZML(J)=0 ODOO 
DO 550 I =2, 1 VI 

ZML(J) =ZML(J) +XF (I, J)*RHOL*VOL(l,J) 

550 CONTINUE 

ZMS = ZMT-ZM L(J) ^ .. 

ZRM = DSQRT(1 . 0D00*RI2*RI2 + ZMS/(RHOS*PI DZZ)) 

ZCHARL(J) = RV-ZRM 
IF (ZCHARL(J) LE 0 ODOO) GO TO 560 
08= DZZ/ZCH ARL(J) 

I = I VI 

IF (RV EQ R01) 1=11 

C NOTE FOR CASES WITHOUT VOID, I— >11 AND T HOT = T SHELL « O.D. 
C 

ZRA(J)=C7*(T(I,J)-TM)*ZCHARL(J)**3 

IF (ZRA(J) LE 0 ODOO) GO TO 560 . 

ZNU(J)=0 42D00*PR**0 012DOO*ZRA(J)**0.25DOO*C18**(— 0.3D00) 

IF (08 LT 10 ODOO) ZNU(J)=0 22D00*(ZRA(J)*PR/(0 2D00+PR)) 

Si **0 28D00*C18**(-0 25D00) 

IF (ZNU(J) LT.l ODOO) GO TO 560 
GO TO 561 

560 ZNU(J) = 1 ODOO 

561 ZKLE(J) = KL*ZNU(J) 

540 CONTINUE 

RETURN 

END 

SUBROUTINE ENERGY(EBALP.TIMEM) 

CHARACTERM PHASE(40,20),NC 

REAL MDOT I KW 1 KS ) KL,KLE,MT I MS,ML,K(40 J 20) 1 KIP(40 1 20) I 

&MT MS VOLS I ML ) VOLL 1 ESHSUM,EFSUM,RA,NU,CHARL I C7 1 C12 ECON(20), 

*TV,TO(40,20) ,T(40,20) ,K,RHO(40,20) ,XF (40,20) 

&YF (40,20) , DR (40) , DZ (20) , Z (20) ,R(40) , EO (40,20) ,E(40,2^,AREA(30), 

&TFO(20) ,TF(20),C1 (40,20) ,C2(40, 20) ,C3(40,20) ,C4(40,20)RCMLO( ) 

ATFOUT,EE(40,20),TIMEPR1,TIMEPR2,TIMEP,TIMEF 1 TFP TFF ^rv.rm, 

&TIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL I EEL(40,20) J 

&EELO(40,20),EEL02(40,20),PR,EML(40,20), 

&MFL.QF, II, IV1,IV2,JJ,N, PHASE, NC, 

&ZKLE (1 0) ,ZNU (10) ,ZRA(1 0) ,ZCHARL(1 0) .ZNUAVG 
QSHELL=Q*C12 
QSALT=QSHELL-QF 
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ESHSUM=QSHELL*DT*(N-1) 

N546 = 3276 DOO/DT 
ESH546=Q1*3276.D00*C12 
N 728 = 4368 DOO/DT 
ESH728=Q2*1092D0O*C12+ESH546 

IFfTIMEM GT 54 7 AND TIMEM LT 72 7) ESHSUM = ESH546+ QSHELL*DT* 

& (N-N546) 

IFfTIMEM GT 72 8) ESHSUM=ESH728+QSHELL*DT*(N-N728) 

ECAN=OOODOO 
DO 835 1 = 1,11 
DO 840 J=1,JJ 
ECAN=ECAN+E(I,J) 

840 CONTINUE 
835 CONTINUE 

EBAL=ESHSUM-EFSUM-(ECAN-EINIT) 

EBALP=100 ‘EBAL/ESHSUM 

RETURN 

END 

SUBROUTINE WALLFRACfWALLFRl ,WALLFR2,WALLFR3) 

CHARACTERM PHASE(40,20),NC 

REAL M DOT, KW, KS, KL, KLE , MT, MS, M L, K(40, 20) , KIP(40,20) , 
8cKIM(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU j NUL,KF 
COMMON DT.DRR.DZZ.DELWO^ELWI.DELWS.Rll.KIP.BBINVflS^.QRADflS), 
&RI2,R01,R02,Q,TFI,MD0T,CPF,KW,CPW I RH0W,K5,CPS,RH0S,FF(18,18), 
&KL,CPL,RHOL,TM,HM,KLE,Z2,Z3,Z4,PI,VOLT,Q1,Q2,VOL(40,20) I RCML(20), 

8tMT, MS, VOLS, ML, VOLL, ESHSUM, EFSUM,RA,NU,CHARL ( C7,C12,ECON(20), 

ATV,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20) J AREAOLD(3), 

*YF(40,20),DR(40),DZ(20),Z(20),R(40),EO(40,20),E(40,20) 1 AREA(30), 

*TFO(20),TF(20),C1(40,20),C2(40,20),C3(40,20),C4{40,20),RCMLO(20) I 

8£TFOUT,EE(40,20),TIMEPR1,TIMEPR2,TIMEP,TIMEF 1 TFP,TFF,RV I RM, 

&TIME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL,EEL(40 I 20), 

&EELO(40,20),EEL02(40,20),PR,EML(40,20), 

&MFL,QF, II, IV1,IV2,JJ,N, PHASE, NC, 
AZKLE(10),ZNU(10),ZRA(10),ZCHARL(10),ZNUAVG 
112 = 11/2 
QW1=0 0D00 
QW2=0 0D00 
QW3=O.ODOO 
QPCM1=0 0D00 
QPCM2=OODOO 
QPCM3=0 0D00 

QW1 = KW*2 ODOO*PI*DELWS*(T(3,l)-T(2,l))/ 

* DLOGfl OODOO*R(3)/R(2)) 

QW1 =QW1 + KW*2 ODOO*PI»DELWS*(T(3,JJ)-T(2,JJ))/ 

L DLOGfl OODOO*R(3)/R(2J) 

QW2=KW*2 0DO0*PI*DELWS*fT(ll2+l,l)-T(ll2,l))/ 

L DLOGfl 00D00*R(II2+1)/R(II2)) 

QW2= QW2+ KW*2 0D00*PI*DELWS*fT(ll2+ 1 , JJ)-T(II2,JJ))/ 
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Sc DLOG(1.00D00*R(H2+l)/R(H2)) 
QW3=KW*2.0D00*PI*DELWS*(T(II,1)-T(II-1,1))/ 

Sc DL0G(1 OODOO*R(ll)/R(ll-l)j 

QW3=QW3+KW*2.0D00*PI*DELWS*(T(II,JJ)-T(II-1,JJ))/ 

Sc DL0G(1 OODOO*R(II)/R(II-1)) 

DO 701 J =2,JJ — 1 

QPCM1=QPCM1 + WP(2 1 J)*2.0DOO*PI*DZZ*(T(3,J)-T(2,J))/ 

Sc DLOG (1 . OODOO* R (3)/R (2)) 

QPCM2=QPCM2+KIP(ll2,J)*2.0D0O*PI*DZZ*fT(ll2+l p J)— T(II2,J))/ 

Sc DL0G(1 00D00*R(II2+1)/R(II2)) 

QPCM3=QPCM3+DABS(ECON(J)/DT)+DABS(QRAD(J+7)) 

701 CONTINUE 

WALLFR1 = 1 00 D00*DABS(QW1)/(DABS(QW1) + DABS(QPCM1)) 

WALLFR2 = 1 00 D00*DABS(QW2)/(DABS(QW2) + DABS(QPCM2)) 

WALLFR3 = 1 00 DOO* DABS (QW3)/(DABS (QW3) + DABS(QPCM3)) 

RETURN 

END 

SUBROUTINE 0UTPUT(TIMEM,N1,N2,WALLFR1 J WALLFR2 I WALLFR3,EBALP I WF) 

CHARACTERM PHASE(40,20),NC 

REAL MDOT.KW.KS.KL.KLE.MT.MS.ML.K^OJ.KIP^O.ZO), 
&KIM(40,20),KJP(40,20),KJM(40,20),MFL,MUF,NU,NUL 1 KF 
COMMON DT.DRR.DZZ.DELWO.DELWI.DELWS.Rll.KIP.BBINVtlMfyQRADllS), 
4iRI2,R01,R02,Q,TFI,MD0T 1 CPF J KW,CPW,RH0W,KS,CPS 1 RH0S,FF(18,18), 
tKL,CPL,RHOLJM 1 HM,KLE,Z2 1 Z3 ) Z4 ) PI,VOLT 1 Ql 1 Q2,VOL(40 l 20),RCML(20) ) 
8 i MT J MS,VOLS 1 ML,VOLL,ESHSUM 1 EFSUM 1 RA,NU 1 CHARL 1 C7 1 C12 > ECON(20), 
A:TV,TO(40,20),T(40,20),K,RHO(40,20),XF(40,20),TFIN(20),AREAOLD(3), 

&YF (40,20) l DR(40),DZ(20) 1 Z(20) J R(40),EO(40 l 20) I E(40 1 20) l AREA(30) 1 
4£TFO(20),TF(20),C1(40,20) 1 C2(40 1 20),C3(40,20) 1 C4(40,20) 1 RCMLO(20) 1 
&TFOUT,EE(40,20),TIMEPR1,TIMEPR2,TIMEP 1 TIMEF,TFP I TFF,RV,RM 1 
8J1ME2,TIME3,TIME4,TF2,TF3,TF4,EINIT,ECAN,EBAL 1 EEL(40 J 20), 

&EELO(40,20),EEL02(40,20),PR,EML(40,20), 

&MFL.QF, II, IV1,IV2,JJ,N, PHASE, NC, 

&ZKLE(10),ZNU(10),ZRA(10),ZCHARL(10),ZNUAVG 

IF(N NE Nl) GO TO 810 

OUTPUT WALL FRACTION VALUES 

WRITE (11, 3400) TIMEM,WALLFR1,WALLFR2,WALLFR3 
JJ2=JJ/2 

QVOIDCON=0 0D00 
QVOIDRAD=0 0D00 
QVOIDTOT=0 0D00 
ZNUAVG = OODOO 
ZRAAVG=0 0D00 
DO 811 J=2,JJ-1 

QVOIDCON=QVOIDCON + (-ECON(J)/DT) 
QVOIDRAD=QVOIDRAD+(-QRAD(J-l)) 

QVOIDTOT =QVOIDCON+QVOIDRAD 
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ZNUAVG = ZNUAVG + ZNU (J) 

ZRAAVG = ZRAAVG + ZRA( J) 

811 CONTINUE 

ZRAAVG=ZRAAVG/(JJ-2) 

ZNUAVG =ZNUAVG/(JJ-2) 

C 

C OUTPUT TIME-DEPENDENT TEMPERATURES AND HEAT FLUXES 
C 

WRITE (8,2500) TIMEM,MFL,T(II,JJ2),WF, 
*TFI,TFOUT,QF,EBAL,EBALP,QVOIDCON,QVOIDRAD,QVOIDTOT 
C 

C OUTPUT FREE CONVECTION VALUES 
C 

IF(NC EQ *ON ’) WRITE (9,2600) TIMEM,CHARL,TV,RA,NU 
IF(NC.EQ 'ON ’) WRITE (17,3750) TIMEM,ZNUAVG,(ZNU(J),J=2,JJ-1), 
ti (ZCHARL(J),J=2,JJ-1) 

IF(NC EQ.’ON ') WRITE (18,3755) TIMEM, ZRAAVG, (ZRA(J) , J =2, J J— 1) 
N1 =N1 + (TIMEPR1/DT) 

GO TO 850 
C 

C OUTPUT PCM PHASE DISTRIBUTIONS 
C 

810 WRITE (13,3000) TIMEM, MFL,WF,NU,(Z(J),J=2,JJ-1) 

WRITE (13,3050) 

DO 800 1=2,11-1 

111=11 + 1-1 

WRITE(13,3100) R(II1),(PHASE(II1,J),J=2,JJ-1),K(II1,5) 

800 CONTINUE 
C 

C OUTPUT CANISTER TEMPERATURE DISTRIBUTIONS 
C 

WRITE (12,3450) TIMEM, MFL,WF,NU,(Z(J),J = 1,JJ) 

WRITE (12,3050) 

DO 855 1 = 1,11 
111=11 + 1-1 

WRITE(12,3500) R(II1),(T(II1 1 J),J = 1,JJ) 

855 CONTINUE 

N2=N2+TIMEPR2/DT 

2500 FORMAT(F6 2,2X,F5 3,1X,F6 1,2X,F6 4,2(2X,F6 1),1X,F7.2, 

& 2X,F7 0,2X,F5 2,2X,3(F8 3,2X)) 

2600 F0RMAT(F6.2,1X,F5 3,2X,F6.1,2X,F9 0,2X,F6.3) 

3000 FORMAT(//’ TIME=’,F6.2,’ MFL=',F6 4,’ WF= , ,F6 4,’ NU=’,F6.3 
&//T25,’ CANISTER PCM PHASE MAP'//’ Z, CM =*.18(1X,F6.4)) 

3050 FORMAT(’R, CM’) 

3100 FORMAT(F6 4,7X,8(A4,3X),D11 4) 

3150 FORMAT(//’ SHELL TEMPERATURES AT TIME = \F6.2,' MIN ’// 
ft' R, CM ’,3X,’ T(R,Z=0) ',3X,’ T(R,Z=L) '/) 

3200 FORMAT(F6 4,6X, F6 1 ,8X,F6 1) 

3250 FORMAT(//’ SHELL TEMPERATURES AT TIME = ’.FM,' MIN ’// 
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ti' TFI=',F6 1,’ TF0UT=’,F6.1// 

V Z, CM , J 3X, , T(R=R0,Z) *,4X 1 ’T(R=RI,Z) ’,5X,' TF(Z) ’/) 

3300 FORMATfFefOX.Fe.l.SX.Fe 1,8X,F6 1) 

3400 FORMAT(F7.3 1 3(3X J F7-3) J 2X,2(l2,2X),6(F6 4,2X)) 

3450 FORMAT (//’ TIME=’,F6 2,’ MFL=’,F6 4,* WF= , ,F6.4 1 ’ NU=’,F6 3 
&//T25/ CANISTER TEMPERATURE MAP’//’ Z, CM =\20(2X,F6.4)) 
3500 FORMAT(F6 4,6X,20(F6.1,2X)) 

3750 FORMATiFe^.FSJWfFS.S.lXj.lX^FS.S.lX)) 

3755 FORMAT(F6.2,2X,F9.0 I 2X 1 8(F9.0 1 1X)) 

850 RETURN 
END 
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A2 . 2 Program Variable Definitions 
VARIABLE 


NAME 

TYPE 

DEFINITION 

AL 

REAL*8 

LIQUID PCM THERMAL DIFFUSIVITY, CM2/SEC 

AREA 

REAL*8 

VOID SURFACE ELEMENT AREA ARRAY, CM2 

AS 

REAL*8 

SOLID PCM THERMAL DIFFUSIVITY, CM2/5EC 

BB 

REAL*8 

VOID 'VIEW FACTOR - EMITTANCE" MATRIX 

BETA 

REAL*8 

VOLUMETRIC COEFFICIENT OF THERMAL EXPANSION, 1/K 

CHARL 

REAL*8 

LIQUID PCM CHARACTERISTIC LENGTH, CM 

CPF 

REAL*8 

HE/XE COOLING FLUID SPECIFIC HEAT, J/G/K 

CPL 

REAL*8 

LIQUID PCM SPECIFIC HEAT, J/G/K 

CPS 

REAL*8 

SOLID PCM SPECIFIC HEAT, J/G/K 

CPW 

REAL*8 

CONTAINMENT WALL SPECIFIC HEAT, J/G/K 

Cl 

REAL*8 

ELEMENT CONDUCTOR COEFFICIENT ARRAY 

C2 

REAL*8 

ELEMENT CONDUCTOR COEFFICIENT ARRAY 

C3 

REAL*8 

ELEMENT CONDUCTOR COEFFICIENT ARRAY 

C4 

REAL*8 

ELEMENT CONDUCTOR COEFFICIENT ARRAY 

C7-C18 

REAL*8 

MISCELLANEOUS CONSTANTS 

DD 

REAL*8 

VOID 'TEMPERATURE 4 DIFFERENCE" MATRIX 

DEL 

REAL*8 

KRONECKER DELTA FUNCTION 

DELWI 

REAL*8 

CANISTER INNER WALL THICKNESS, CM 

DELWO 

REAL*8 

CANISTER OUTER WALL THICKNESS, CM 

DELWS 

REAL*8 

CANISTER SIDE WALL THICKNESS, CM 

DR 

REAL*8 

RADIAL GRID SIZE ARRAY, CM 

DT 

REAL*8 

TIME STEP, SEC 

DUM 

REAL*8 

DUMMY VARIABLE 

DZ 

REAL*8 

AXIAL GRID SIZE ARRAY, CM 

E 

REAL*8 

N + l TIME STEP ELEMENT ENTHALPY, J 

EBAL 

REAL*8 

GLOBAL ENERGY BALANCE, J 

ECAN 

REAL*8 

CANISTER ENERGY CONTENT, J 

ECON 

REAL*8 

VOID ELEMENT CONDUCTION ENERGY ARRAY, J 

EE 

REAL*8 

ELEMENT SPECIFIC ENTHALPY, J/G 

EFSUM 

REAL*8 

SUM OF ENERGY TRANSFER TO COOLING FLUID, J 

EINIT 

REAL*8 

INITIAL CANISTER ENERGY CONTENT, J 

EO 

REAL* 8 

N TIME STEP ELEMENT ENTHALPY, J 

EPS 

REAL*8 

VOID SURFACE ELEMENT EMITTANCE ARRAY 

ES 

REAL*8 

PCM EMITTANCE 

ESHSUM 

REAL*8 

ABSORBED CANISTER ENERGY, J 

EW 

REAL*8 

CONTAINMENT WALL EMITTANCE 

FF 

REAL*8 

VOID SURFACE ELEMENT VIEW FACTOR ARRAY 

FFSUM 

REAL*8 

SUM OF VOID SURFACE ELEMENT VIEW FACTORS 

G 

REAL*8 

GRAVITY ACCELERATION, CM/SEC2 

GAMMAO 

REAL*8 

1-RHOL/RHOS 

H 

REAL*8 

COOLING FLUID FILM COEFFICIENT, W/CM2/K 

HM 

REAL*8 

PCM HEAT OF FUSION, J/G 

1 

INT*4 

RADIAL DO LOOP INDEX 

II 

INT*4 

TOTAL NUMBER OF RADIAL GRIDS 
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VARIABLE 


NAME 

TYPE 

IV1 

INT*4 

IV2 

INTM 

J 

INTM 

JJ 

INTM 

K 

REALM 

KIM 

REALM 

KIP 

REALM 

KJM 

REALM 

KJP 

REALM 

KF 

REALM 

KL 

REALM 

KLE 

REALM 

KS 

REALM 

KW 

REALM 

MDOT 

REALM 

MFL 

REALM 

ML 

REALM 

MS 

REALM 

MT 

REALM 

MUF 

REALM 

N 

INTM 

NC 

CHAR 

NN 

INTM 

NRS 

INTM 

NU 

REALM 

NUL 

REALM 

N1,N2 

INTM 

PHASE 

CHAR 

PI 

REALM 

PR 

REALM 

PRF 

REALM 

Q 

REALM 

QF 

REALM 

QPCM 

REALM 

QRAD 

REALM 

QSALT 

REALM 

QSHELL 

REALM 

QVOIDCON 

REALM 

QVOIDRAD 

REALM 

QVOIDTOT 

REALM 

QW 

REALM 

R 

REALM 

RA 

REALM 

RE 

REALM 

RHO 

REALM 

RHOF 

REALM 


DEFINITION 


RADIAL PCM GRIDS ADJACENT TO VOID 

RADIAL VOID GRIDS ADJACENT TO PCM 

AXIAL DO LOOP INDEX 

TOTAL NUMBER OF AXIAL GRIDS 

ELEMENT THERMAL CONDUCTIVITY ARRAY, W/CM/K 

"1-1" ELEMENT-TO-ELEMENT CONDUCTOR, W/CM/K 

"1 + 1" ELEMENT-TO-ELEMENT CONDUCTOR, W/CM/K 

"J-l" ELEMENT-TO-ELEMENT CONDUCTOR, W/CM/K 

"J+l" ELEMENT-TO-ELEMENT CONDUCTOR, W/CM/K 

COOLING FLUID THERMAL CONDUCTIVITY, W/CM/K 

LIQUID PCM THERMAL CONDUCTIVITY, W/CM/K 

ENHANCED LIQUID PCM THERMAL CONDUCTIVITY, W/CM/K 

SOLID PCM THERMAL CONDUCTIVITY, W/CM/K 

CONTAINMENT WALL THERMAL CONDUCTIVITY, W/CM/K 

COOLING FLUID MASS FLOW RATE, G/SEC 

MASS FRACTION LIQUID PCM 

TOTAL LIQUID PCM MASS, G 

TOTAL SOLID PCM MASS, G 

TOTAL PCM MASS, G 

COOLING FLUID VISCOSITY, G/SEC/CM 

TIME STEP DO LOOP INDEX 

NATURAL CONVECTION ON/OFF FLAG 

TOTAL NUMBER OF TIME STEPS 

NUMBER OF VOID SURFACE RADIATING ELEMENTS 

LIQUID PCM NUSSELT NUMBER 

LIQUID PCM KINEMATIC VISCOSITY, CM2/SEC 

UPDATE AND/OR PRINTED OUTPUT TIME STEPS 

PCM ELEMENT PHASE, LIQUID, MUSHY, SOLID, OR VOID 

PI CONSTANT 

LIQUID PCM PRANDTL NUMBER 

COOLING FLUID PRANDTL NUMBER 

OUTER WALL ABSORBED HEAT FLUX, W/CM2 

HEAT TO COOLING FLUID, W 

RADIAL HEAT TRANSFER IN PCM, W 

VOID SURFACE ELEMENT NET HEAT LOSS ARRAY, W 

NET HEAT TO CANISTER, W 

OUTER WALL ABSORBED HEAT, W 

VOID CONDUCTION HEAT TRANSFER, W 

VOID RADIATION HEAT TRANSFER, W 

VOID TOTAL HEAT TRANSFER, W 

RADIAL HEAT TRANSFER IN CANISTER SIDE WALLS, W 

RADIAL COORDINATE, CM 

LIQUID PCM RAYLEIGH NUMBER 

COOLING FLUID REYNOLDS NUMBER 

ELEMENT DENSITY ARRAY, G/CM3 

COOLING FLUID DENSITY, G/CM3 
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VARIABLE 


NAME 

TYPE 

DEFINITION 

RHOL 

REAL*8 

LIQUID PCM DENSITY, G/CM3 

RHOS 

REAL*8 

SOLID PCM DENSITY, G/CM3 

RHOW 

REAL*8 

CONTAINMENT WALL DENSITY, G/CM3 

Rll 

REAL*8 

COOLING FLUID TUBE INNER RADIUS, CM 

RI2 

REAL*8 

Rll + DELWI, CM 

RM 

REAL*8 

AVERAGE PCM SOLID-LIQUID INTERFACE RADIUS, CM 

ROl 

REAL*8 

R02-DELW0, CM 

R02 

REAL*8 

CANISTER OUTER RADIUS, CM 

RV 

REAL*8 

PCM-VOID INTERFACE RADIUS, CM 

SIGMA 

REAL*8 

STEFAN-BOLTZMANN CONSTANT, W/CM2/K4 

T 

REAL*8 

N + l TIME STEP ELEMENT TEMPERATURE, K 

TF 

REAL*8 

COOLING FLUID TEMPERATURE ARRAY, K 

TFI 

REAL*8 

COOLING FLUID INLET TEMPERATURE, K 

TFOUT 

REAL*8 

COOLING FLUID OUTLET TEMPERATURE, K 

TIMEM 

REAL*8 

SIMULATION TIME, MIN 

TM 

REAL*8 

PCM MELTING POINT, K 

TO 

REAL*8 

N TIME STEP ELEMENT TEMPERATURE, K 

TSHELL 

REAL*8 

OUTER CANISTER WALL AVERAGE TEMPERATURE, K 

TTUBE 

REAL*8 

AVERAGE COOLING FLUID TUBE TEMPERATURE, K 

7V 

REALM 

INNER VOID SURFACE TEMPERATURE, K 

U 

REALM 

OVERALL HEAT TRANSFER COEFFICIENT, W/CM2/K 

V 

REALM 

COOLING FLUID VELOCITY, CM/SEC 

VOL 

REALM 

ELEMENT VOLUME ARRAY, CM3 

VOLL 

REALM 

TOTAL LIQUID PCM VOLUME, CM3 

VOLS 

REALM 

TOTAL SOLID PCM VOLUME, CM3 

VOLT 

REALM 

TOTAL VOLUME, CM3 

VOLV 

REALM 

TOTAL VOID VOLUME, CM3 

WF 

REALM 

VOID VOLUME FRACTION 

WALLFR1 

REALM 

RADIAL HEAT TRANSFER WALL FRACTION AT INNER RADIUS 

WALLFR2 

REALM 

RADIAL HEAT TRANSFER WALL FRACTION AT MEAN RADIUS 

WALLFR3 

REALM 

RADIAL HEAT TRANSFER WALL FRACTION AT OUTER RADIUS 

XF 

REALM 

ELEMENT LIQUID PCM MASS FRACTION 

YF 

REALM 

ELEMENT LIQUID PCM VOLUME FRACTION 

Z 

REALM 

AXIAL COORDINATE, CM 

ZCHARL 

REALM 

AXIAL DEPENDENT CHARACTERISTIC LENGTH ARRAY, CM 

ZKLE 

REALM 

AXIAL DEPENDENT ENHANCED CONDUCTIVITY ARRAY, W/CM/K 

ZNU 

REALM 

AXIAL DEPENDENT NUSSELT NUMBER ARRAY 

ZRA 

REALM 

AXIAL DEPENDENT RAYLEIGH NUMBER ARRAY, CM 

ZRM 

REALM 

AXIAL DEPENDENT SOLID-LIQUID INTERFACE POSITION, CM 
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Appendix A3 . Video Ani mations 

A video tape which animates the transient, numerical 
results from two-dimensional canister analyses was created 
at the NASA Lewis Research Center (LeRC) advanced Graphics 
Vis ualization Lab oratory (G-VIS Lab). The animation 
visually depicts canister temperatures, temperature 
gradients, and PCM phase distributions through the combined 
use of color fringe and isotherm contour plotting 
techniques. A data set containing temperature predictions 
for a single 91 minute TES charge-discharge cycle at 6- 
second intervals is read in and animated by the LeRC- 
developed program SVP (Scientific Visualization Program). 

SVP is run on the LeRC VM mainframe system through a Silicon 
Graphics IRIS4D/120 workstation. The graphical output from 
SVP is transferred to an Abekas A60 Digital Video Disk 
Recorder with which animation loops and segments are defined 
and displayed at various speeds. In conjunction with the 
A60, an Abekas A34 Solo unit is used for editing and video 
special effects before the final video sequences are 
transferred to 1 inch video tape, 3/4 inch (Umatic) or VHS 
videocassette for presentation. This procedure is repeated 
for each of three cases of two-dimensional canister 
analyses: without a void model, with a void model, and with 

void and free convection models. 

The animated numerical results are displayed in real 
time at 100:1 or 300:1 time compression ratios. For 
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example, temperature predictions from a 91 minute cycle 
displayed at a 300:1 time compression ratio would have a 
viewing time of 91/300 minutes or 18.2 seconds. This kind 
of visualization provides an effective method to confirm the 
accuracy of numerical predictions, more thoroughly interpret 
spatial and temporal relationships, and observe complex 
phenomena which can not be observed in experiments. From a 
practical standpoint, visualization of numerical predictions 
reduces the task of reviewing a 2-inch thick computer output 
listing containing 182,000 tabular temperature predictions 
(10 by 20 finite-difference elements x 91 minute cycle x 10 
up-dates per minute = 182,000 predictions) to viewing an 
18.2 second video tape presentation. 

The 12-minute, VHS video tape animating two-dimensional 
canister temperature predictions is available for viewing by 
contacting the author or the Cleveland State University 
Department of Mechanical Engineering (Dr. Mounir B. 

Ibrahim) . 
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